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ABSTRACT 

In  this  thesis,  an  adaptive  lattice  algorithm  is  derived  for  an  ARM  A  digital  lattice 
filter,  whose  parameters  are  estimated  using  a  generalized  Mullis-Roberts  criterion  for 
parameter  estimation.  Design  of  the  ARMA  lattice  filter  based  on  this  generalized  cri- 
terion is  studied  as  is  the  accuracy  of  the  parameter  estimation  algorithm  used  in  its  de- 
sign. Application  of  the  derived  lattice  algorithm  to  system  identification  modeling  is 
demonstrated  through  computer  simulation  of  various  system  identification  problems. 
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I.     INTRODUCTION 

Digital  signal  processing  is  a  field  which  is  rapidly  expanding  due  to  advances  in 
modern  technology.  Essential  to  this  field  are  digital  filters.  Modeling  these  filters  con- 
stitutes much  of  the  effort  involved  in  digital  signal  processing.  The  filters  provide  a 
transfer  function  which  describes  the  relationship  between  filter  input  and  output. 
Autoregressive  (AR),  moving  average  (MA)  and  combination  autoregressive  moving 
average  (ARMA)  models  are  widely  used  to  represent  the  transfer  function  of  a  digital 
filter.  A  filter  transfer  function  is  commonly  described  in  direct  form.  This  form  is  a  ratio 
of  two  polynomials,  usually  of  the  form, 

B(z)  b^l^z'1  +  -  +  blz~c 

A^)  1  +  al  z      +  •••  +  a5z 

The  above  equation  describes  an  ARMA  model  of  order  (s,t)  where  s  is  the  order  of  the 
denominator  and  t  the  order  of  the  numerator.  The  as  parameters  form  the 
autoregressive  portion  of  the  ARMA  model.  The  b,  parameters  form  the  moving  average 
portion  of  the  ARMA  model.  If  all  the  autoregressive  parameters  are  zero,  then  the  filter 
transfer  function  H(z)  is  strictly  a  moving  average  process  of  order  t.  If  all  the  moving 
average  parameters  are  zero  except  for  b0  equal  to  one,  then  the  filter  transfer  function 
is  strictly  an  autoregressive  process  of  order  s. 

A.     OBJECTIVES  OF  THE  THESIS 

Fundamental  to  the  design  of  digital  filters  is  estimation  of  AR,  MA  or  ARMA 
parameters.  Accurate  and  efficient  parameter  estimation  has  been  the  subject  in  much 
of  the  related  digital  signal  processing  literature  [Refs.  1  ,2.3].  The  first  objective  of  this 
thesis  is  to  confirm  the  proposed  ARMA  parameter  estimation  algorithm  of  [Ref.  4:  pp. 
619-621],  which  leads  to  the  design  of  a  new  ARMA  digital  lattice  filter.  The  proposed 
algorithm  is  a  generalization  of  the  Mullis-Roberts  criterion  for  parameter  estimation 
known  as  the  modified  least  squares  problem  [Ref.  5:  pp.  227-228].  The  algorithm  uses 
two  recursive  formula  to  estimate  the  parameters.  One  is  an  AR  recursive  formula  which 
estimates  ARMA  parametes  as  the  AR  order  is  increased  by  one.  The  other  is  an  MA 
recursive  formula  which  estimates  ARMA  parameters  as  the  MA  order  is  increased  by 
one.  This  algorithm  is  unique  in  that  it  allows  for  the  design  of  an  ARMA  model  with 


arbitrary  AR  and  MA  orders  with  no  dependency  of  an  AR  model  on  an  MA  model  or 
vice  versa.  The  ARMA  digital  filter  designed  from  the  proposed  ARMA  parameter  es- 
timation algorithm  is  in  the  form  of  a  lattice  structure.  Lattice  realizations  of  filters  are 
widely  used  and  provide  excellent  analysis  of  prediction  errors  [Refs.  6  :  pp.  165-168,7]. 
Gray  and  Markel  developed  an  algorithm  which  produces  lattice  realizations  of  filters 
from  the  direct  form  [Ref.  8]. 

The  second  objective  of  the  thesis  is  to  make  the  proposed  ARMA  digital  lattice 
filter  of  [Ref.  4:  p.  662]  adaptive.  An  adaptive  lattice  filter  is  one  in  which  the  lattice 
coefficients  are  automatically  adjusted  by  an  adaptive  algorithm  to  yield  the  optimum 
filter  design.  The  adaptive  lattice  algorithm  derived  in  this  thesis  is  based  on  the  widely 
used  least  mean  square  (LMS)  algorithm.  Adaptive  filters  have  many  applications  [Ref. 
9:  pp.  7-31]  including. 

1.  System  identification. 

2.  Digital  representation  of  speech. 

3.  Adaptive  auotoregressive  spectrum  analysis. 

4.  Adaptive  detection  of  a  signal  in  noise  of  unknown  statistics. 

5.  Echo  cancellation. 

6.  Adaptive  line  enhancment. 

7.  Adaptive  beamforming. 

The  need  for  an  adaptive  filter  is  made  apparent  by  considering  a  filter  of  fixed  design 
which  is  optimized  for  given  input  conditions.  In  practice,  the  complete  range  of  input 
conditions  may  not  be  known  or  could  change  from  time  to  time.  A  filter  of  fixed  design 
would  not  produce  optimum  results  under  these  conditions.  An  adaptive  filter,  which 
yields  optimum  results  given  changing  input  conditions,  will  give  superior  performance 
to  one  of  fixed  design. 

The  last  objective  is  to  analyze  convergence  properties  of  the  derived  adaptive  lattice 
algorithm.  This  is  accomplished  by  computer  implementation  of  the  adaptive  algorithm. 
The  output  of  a  known  transfer  function  is  compared  to  the  output  of  the  adaptive 
lattice  filter  given  a  common  input.  Plots  of  the  error  between  the  two  outputs  and 
lattice  coefficient  convergence  are  obtained. 

B.     ORGANIZATION  OF  THE  THESIS 

Chapter  II  is  designed  to  present  the  ARMA  parameter  estimation  algorithm  and 
ARMA  digital  lattice  filter  proposed  in  [Ref.  4:  pp.  617-628].  Computer  simulation  of  the 


algorithm  was  performed  and  results  are  shown.  A  brief  review  of  the  Mullis-Roberts 
criterion  is  provided  to  establish  a  reference  for  expanding  this  criterion  to  the  proposed 
AR.MA  parameter  estimation  algorithm.  An  adaptive  lattice  algorithm  is  derived  in 
Chapter  III  which  makes  the  proposed  ARMA  digital  lattice  filter  adaptive.  The  adap- 
tive lattice  algorithm  is  efficient  and  accurate.  Chapter  IV  contains  experimental  results 
which  show  convergence  aspects  of  the  adaptive  lattice  algorithm  when  applied  to 
ARMA  lattice  filters.  Conclusions  about  the  proposed  ARMA  parameter  estimation 
algorithm  as  well  as  the  derived  adaptive  algorithm  are  discussed. 


II.     ARMA  DIGITAL  LATTICE  FILTER 

In  this  chapter,  we  will  review  the  Mullis-Roberts  criterion  for  solving  linear  ap- 
proximation problems  and  introduce  analysis  equations  of  the  ARMA  digital  lattice  fil- 
ter. The  criterion  used  in  the  formulation  of  the  ARMA  digital  lattice  filter  is  a 
generalized  form  of  the  Mullis-Roberts  criterion  [Ref.  5:  pp.  227-22S],  which  has  been 
given  as  a  modified  least  mean  square  problem  for  ARMA  parameter  estimation. 

A.     MULLIS-ROBERTS  CRITERION 

The  Mullis-Roberts  criterion  evolved  from  considering  second  order  statistics  in 
conjuction  with  first  order  information  about  a  process  to  obtain  filter  approximations. 
Consider  the  bounded  impulse  response  sequence  h  =  {/?0,  hu  ...  }  containing  first  order 
information  about  the  filter  h  having  a  frequency  response  function, 

oo 

H{^)  =  Yjhke-JMk  (2.1) 

Let  {  u,  }  be  a  zero-mean,  unit-variance,  white-noise  sequence  and  {  y,  }  be  the  output 
process  corresponding  to  the  input  {  ut  }  ,  then  we  have  the  following  convolutional  re- 
lationship given  by, 

oo 

A'=0 

Second  order  information  about  the  filter  h  is  obtained  from  the  autocorrelation  se- 
quence {  rk  }  given  as 

oo 

>*  =  E{ytyt+^  =  2_jhihk+t  (2-3) 

(=0 

From  equations  (2.1)  and  (2.2),  the  second  order  interpolation  problem  is  to  find  a 
lowest  order  recursive  filter  which  matches  the  data  {  /i0, ...  ,  /*„,.  r0l ...  ,  rm  }  ,  where  h,  re- 
presents the  first  order  information  and  r]  the  second  order  statistics. 


Let  us  now  consider  the  case  where  only  first  order  information  about  a  process  is 
known.  That  is,  given  an  impulse  response  sequence  {  fy,,  hu  ...  }  ,  we  want  to  find  a  re- 
cursive filter  of  the  form, 


(2.4) 
1  +  a}z      +  ••■  +  anz 


k=o 


which  approximates  H(z)  and  therefore  the  impulse  response  sequence  {  h0.  hu  ...  }  .  We 
also  desire  the  frequency  response  H(ejw)  to  approximate  the  desired  response  H(eJm). 
Suppose  that  //(ew)  is  chosen  such  that  it  minimizes  the  integral  squared  error, 


~)tt 


|//(0  -  H(eiM)\2^     =     \\h-hf  (2.5) 


Using  the  Parseval  relation,  we  can  obtain  an  alternative  definition  of  the  approximation 
error. 


k=Q 

If  the  filters  H(z)  and  H{z)  are  driven  by  the  same  white  noise  source,  equation  (2.6) 
describes  the  output  error  between  the  filters  which  we  write  as, 

\\h-hf     =     E(yt-y£     =     E(ef  (2.7) 

where  y,  and  y,  are  the  outputs  of  the  respective  filters  when  driven  by  the  same  white 
noise  source  as  in  (2.2).  Minimizing  (2.7)  is  a  nonlinear  programming  problem  requiring 
the  entire  impulse  response  sequence.  As  a  result,  computational  efforts  for  obtaining  a 
solution  are  inefficient.  A  modification  to  the  problem  was  introduced  [Ref.  5:  pp. 
227-228]  which  considered  a  cost  function  that  is  quadratic  in  the  coefficients  of  the  re- 
cursive filter  given  by  (2.4).  The  modification  is  described  as  follows.  Let 

Q(z)  =  zN(qQ  +  qxz~X  +  -  +  qmz~m)  (2.8) 


ana 


A(=)  =  z\\  +a]Z   '  +  -  +  anz   n) 


(2.9) 


be  the  numerator  and  denominator  polynomials,  respectively,  in  (2.4),  where  A'=  max 
(nun).  The  task  now  is  to  find  coefficients  which  minimize  the  quadratic  form. 


W  =  TT 


\nCJ<")A{c>w)-0{e>M)\7du 


(2.10j 


This  is  a  standard  quadratic  minimization  problem  whose  integral  can  be  expressed  in 
terms     of     the     coefficients     of     polynomials     A(z)     and     Q(z)     and     the     data 

{f\y hm,  r0.  ...  ,  )■„}.  relating  to  the  filter  HC).  This  problem  is  shown  in  Figure  1  and 

equation  (2.10)  is  known  as  the  Mullis-Roberts  criterion. 


u1 

t 

et 

H(z) 

A(z) 

■4> 

► 

Q(z) 

Figure   1.      Modified  least  squares  problem 


B.     GENERALIZED  MULLIS-ROBERTS  CRITERION 

In  order  to  define  the  new  criterion  used  for  ARM  A  parameter  estimation  we  con- 
sider the  following  transfer  function  with  input  sequence  {  x{k),  k  =  1,2, ...  }  and  output 
sequence  {y(k),  k  =  1,2,  ...  }  written  as, 


H(z)  = 


Hy(z) 


(2.11) 


x(k) 

u(k) 

y(k) 

1/H       (2) 
X 

Hy(z) 

Figure  2.      Equivalent  input/output  model 

where  Ht{z)  and  H,{z)  are  reference  polynomials  which  we  desire  to  model.  An  equivalent 
model  is  shown  in  Figure  2  where  u(k)  is  an  intermediate  signal,  and  the  realization  is 
similar  to  that  of  the  direct  form  realization  II  [Ref.  10:  p.  151].  Let  the  intermediate 
sequence  {  n{k),  k  =  1.2, ...  }  be  a  zero-mean  white  gaussian  process.  The  model  of  Fig- 
ure 2  can  then  be  transformed  into  the  model  of  Figure  3  with  u(k)  as  the  common  input 
to  both  transfer  functions  HJz)  and  H,{z). 


u   (k) 

y(k) 

Hy(z) 

X(k) 

Hx(z) 

Figure  3.      Transformed  model 


Note  that  we  earlier  defined  transfer  functions  Hv{z)  and  H„{z)  as  polynomials  of  the 
reference  model.  For  each  transfer  function  an  estimation  polynomial  is  defined  such 


that 


A(z)  =  1  -f  a\z      +  ■■■  +  asz  '      estimates    IIx(z) 


-l 


B(kz)  =  b0  +  blz      +  -  +  btz         estimates    Hv{z) 


(2.12) 
(2.13) 


where  a,  and  b„  (i  =  1.2 s)  and  (/'  =  1,2 /)  ,  are  the  AR  and  MA  parameters,  re- 
spectively, of  the  combined  ARMA  model  formed  by  A(z)  and  B^z).  This  refined  least 
squares  problem  is  shown  in  Figure  4  and  is  a  generalized  form  of  the  Mullis-Roberts 
criterion.  If  the  reference  model  polynomial  //,(r)  in  Figure  4  is  equal  to  unity,  we  obtain 
the  Mullis-Roberts  criterion  shown  in  Figure  1.  Therefore,  by  including  reference 
polynomial  Hx{z),  the  new  criterion  for  ARMA  parameter  estimation  becomes, 


a, 
E    = — — 
■'      2ni 


\Hv(z)A(z)-Hx(z)B(z)\ 


2  dz 


(2.14) 


Figure  4.      Refined  least  squares  problem 


Minimizing  £ii(  is  accomplished  by  calculating  the  coefficients  of  A(z)  and  B(z)  which 
minimize  e(k):  of  Figure  4.  Another  form  of  eq  (2.14)  is  obtained  by  applying  Parseval's 
theorem  and  is  expressed  as, 


ESJ  =  E[(A(z)y(k)-B(z)x(k))2l 


(2.15) 


which  is  obvious  from  Figure  4.    The  coefficients  a„  (i=  1 s),  and  b},  (/' =  1, ...  ,  t), 

which  minimize  £jr  can  then  be  calculated  using  the  normal  equations  for  ARMA 


parameter  estimation.  In  order  to  obtain  the  normal  equations  for  the  problem  in  (2.15), 
let  us  define  the  following: 

av=  [a,  ...  aj    and    bst=  [ft,  ...  bt~\  (2.16) 

are  the  vectors  of  AR  and  MA  parameters,  respectively, 

\t{k)  =  [  j(A)  ...y(k  -  s)    -x(k) ...  -x(k  -  i)  ]  (2.17) 

is  the  data  vector  consisting  of  both  input  and  output  data  elements  and 

R5i,=  £[hv(*)Xr(*)]  (2-18) 

is  the  data  autocorrelation  matrix.   The  criterion  is  to  minimize  the  mean  squared  error 

ESJ  =  E  [  e\k)  ]  =  E[  [  [  1  aSJ  b0  bSJ  ]  h£  ]2] 

rr  r  r  121  ^2'19) 

=  £[[>-(/c)+  [ai>f  60  b.Jh^J-J 

Now  following  the  standard  calculus  of  variation  optimization  procedure  for  minimizing 
EStt,  [Ref.  11:  pp.  100-110],  yields  the  normal  equations  in  the  matrix  form, 

[  1  av   b0  bSJ  ]  R,,r=  [  mm  ESJ  0   0  0]  (2.20) 

It  is  interesting  to  note  that  if  ESJ  in  equation  (2.14)  is  zero,  then  we  have  the  following 
equality, 

Hy(z)        B(z) 

so  the  estimate  for  the  total  reference  model  H(z)  is  the  ratio  of  B(z),  the  MA  part  and 
A(z).  the  AR  part  of  an  estimated  ARMA  model. 

C.     ARMA  PARAMETER  ESTIMATION 

Now  that  the  criterion  for  ARMA  parameter  estimation  has  been  established,  the 
solution  method  to  estimate  the  ARMA  model  parameters  to  minimize  equation  (2.14) 
is  considered.  Let  x(k)  and  y(k)  be  the  input  and  output  signals,  respectively,  of  the  es- 
timated ARMA  model.  Using  a  difference  equation  representation,  this  process  is  de- 
scribed bv 


J  I 

y(k)  =  -Vaj  v(A  -j)  +  Y  bt  x(k  -  i)  (2.22) 

;'=1  (=0 

For  these  input  and  output  signals  we  define  four  estimation,  or  prediction,  models  as 
follows.  The  forward  estimation  signal  for  x(k)  is  defined  as, 


i  j 

)j(k)  =  -Yb*  x(k  -  i)  +  Y  af  y(k  -j)  (2.23) 


7=1 


where  b',  and  a*}  are  the  corresponding  estimation  parameters.  The  forward  esitmation 
signal  for  >•(/%)  is  similarly  defined  as, 


yj{k)  =  -Ya}  y(k  -j)  +  £  bf  x(k  -  i)  (2.24) 


7=1  i=l 

The  backward  estimation  errors  for  x(k  —  /)  and  y(k  —  s)  are  then  given  by, 


{-]  r-l 

xb(k  -  /)  =  -  y  bf  x(k  -  i)  +  y  af  y(k  -j)  (2.25) 


i=o  y=o 


5-1  I-\ 


yb(k  -s)  =  ~Yjaf  y(k  -J)  +  Y  bf  x(k  -  i)  (2.26) 

;=o  ;=0 

where  the  superscripts  g  and  d  indicate  the  backward  estimation  parameters  for  x  and 
y.  respectively. 

From  these  estimation  models,  we  can  now  obtain  the  four  prediction  errors  at  any 
given  time  k.  These  errors  are  expressed  as  differences  between  the  predicted  value  and 
actual  value  of  an  input  or  output  signal,  namely. 

vt1(k)  =  -x(k)  +  xJ[k)  (2.27) 

<,(!<)  =y{k)-yj(k)  (2.28) 

gs,(k)  =  x(k-t)-xb{k-t)  (2.29) 
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djk)=y(k-s)-yb(k-s) 


(2.30) 


We  now  use  the  vector  notation  to  simplify  the  expressions  for  prediction  errors.  In  the 
following,  the  forward  error  elements  corresponding  to  both  x  and  y  form  a  vector 
\Jk),  given  by 


V(*)  =  [    -<,(*)    <,(*)]  =  M*)<r 


and  the  backward  error  vectors  are  given  by 


Ss,t(k)  =  -  M*) G;.' 


ds>l{k)  =  hSil(k)  Djtl 


(2.31) 

(2.32) 
(2.33) 


where  Cst  and  GJf  and  Dv  are  the  forward  estimation  parameter  matrix  and  backward 
estimation  parameter  vectors,  respectively,  defined  as 


Cw  = 


0    a 


-*        i      ;  x 


b: 


1    af    ...    c?s      0    £}    ...    b) 


Ger=W     -     «?-i    0      &£     ...     fcjL,     1] 


DV=[« 


a._i     1       On 


°U  o] 


(2.34) 


(2.3; 


(2.36) 


It  can  be  shown  that  the  prediction  errors  satisfy  some  orthogonal  conditions.  These 
conditions  are  similar  to  those  found  in  AR  modeling  problems  [Ref.  6  :  pp.  116-121]. 
We  now  list  the  orthogonality  conditions  for  the  ARMA  formulation  in  discussion 
without  proof  as  the  following: 

E[  «£(*)    v(k-j)  ]=0        £[  i~»  y(k-j)  ]=0 
£[  <f(£)   x(k-i)  ]  =  0       £[  if/A)   *(*-/)  ]=0 


^U,(A--1)  j-(A--7)]=0 
£U,(A--1)  a-(A--/)]=0 
EU/k-l)  y(k-j)j=Q 


(3.37) 
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E\_dSit{k-  1)   x{k-i)1  =  Q 

where  /=  1,2, ...  ,  t  and y  =  1,2, ...  ,  s. 

In  Section  B  we  have  obtained  a  set  of  normal  equations  in  terms  of  RJf  and  the 
ARM  A  estimation  parameters.  In  this  section,  we  have  defined  four  sets  of  forward  and 
backward  estimation  parameters  and  established  some  orthogonal  relationships.  In  what 
follows,  we  derive  a  set  of  equations  which  relates  the  coefficients  of  the  estimated 
ARMA  model  with  those  of  the  forward  estimation  parameter  matrix  Cut.  Consider  the 
expected  value  of  the  forward  prediction  error,  vsj(k)  and  the  data  hst(k).  Since  the  pre- 
diction error  is  orthogonal  to  all  past  samples  of  data  y(k  —J),  x(k  —  i)  but  not  to  y(k) 
or  x(k)  as  listed  in  (2.37),  the  result  is  a  matrix  which  is  defined  as 


E[V(A-)rhv(£)]  = 


si 


0    £,    0 


cS    0    cA    0 


<>3 


(2.3S) 


where 


f,  =  -E  [  v*t(k)y(k)  ]  =    -£  [  v,»  r5»  ]  =  Eff 


(2.39) 


is  the  crosscorrelation  between  the  forward  prediction  errors  of  x  and  y  at  lag  zero. 

c2  =  E  [  £(*)  *(*)  ]  =  E [  KM?  ]  =  Kt  (2-40) 


is  the  forward  prediction  error  power  of  x(k) 

=  E[Sjk)y(k)]=  E  [  (v-;»)2  ]  =  Eit 


«*3 


(2.41 


is  the  forward  prediction  error  power  of y{k)  and 


c4  =  -E  [  i&A)  x(k)  ] £  [  iftA)  <,(A)  ]  =  £5V 


(2.42) 


is  the  crosscorrelation  between  the  forward  prediction  errors  of  v  and  x  at  lag  zero.  In 
another  interpretation,  the  left  hand  side  of  equation  (2.38)  can  be  written  as, 


£  [  V(*)rM*)  ]  =  E  [  Cw  h5>,(A-)7'hv(A)  ]  =  Cw  RS>1 


(2-43) 


and  we  have 
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*-v     V 


j7x  y 


0     Es,    0 


.44) 


from  equation  (2.38). 
A  similar  approach  for  both  backward  prediction  errors  yields  the  following 


Gv 


Rw 


0    Eff   0    ^i 
o    4    0   ^/ 


(2.45) 


In  order  to  express  the  coefficients  of  the  ARM  A  estimation  model  in  terms  of  the  co- 
efficients of  the  parameter  estimation  matrix  Cst  and  parameter  vectors  Gs ,  and  Dst  , 
consider  the  combination  of  the  normal  equation  (2.20)  and  the  parameter  estimation 
matrix  and  vectors.  From  equation  (2.44)  the  normal  equation  for  ARMA  parameter 
estimation  mav  be  written  as 


0  a.,,     1     *>SI 

1  aj,   0    bj. 


R,,r 


Exy    0    Ex     0 

-v 


H,    0    EtJ   0 


(2.46) 


Combining  equation  (2.46)  with  (2.20)  we  obtain 


1    *s,:   bo    ht 

0  a£,    1    b£, 

1  av    °    blr 


Rv  = 


£?? 


£J(min  0  0  0 
£?/  0  £J  0 
£*,       0    £^v    0 


(2.47) 


In  equation  (2.47)  multiply  row  2  by  -—-  then  subtract  row  2  from  row  3  and  equate 
the  result  to  row  1.  We  then  obtain  the  following  relations 


7xy 


as,t  —  a"j,r       aj\; 


E  J 


£* 

•^j,/- 


(2.4S) 


h  = 


r*y 

^5,1 


(2.49) 


Kt =  Kt  ~  h*,t 


xy 


E  J 
~Er 


(2.50) 


and, 
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£  ,  min  =  Eys,  - 


(57)2 


(2.51) 


Calculation  of  the  ARMA  estimation  model  coefficients  requires  knowledge  of  the  esti- 
mation parameters  a;r,  a] r,  b;f,  b£„  60  and  the  values  £;„  £>r.  E*j  .  Recursive  formulas 
calculate  the  estimation  parameters  as  the  AR  or  MA  order  of  the  estimation  model  in- 
creases by  one.  Given  the  parameters  of  the  ARMA  estimation  model  of  order  (s  ,  t) 
these  formulas  calculate  the  (s  +  1  ,  /)  or  (s ,  t  +  1)  ARMA  parameters.  For  a  compre- 
hensive derivation  of  these  recursive  formulae  see  [Ref.  4:  pp.  619-621].  The  AR-type 
recursive  formula  for  the  forward  parameter  estimation  matrix  and  backward  estimation 
parameter  vectors  is 


Gs+\,i=  [  Gs,t  +  W2D;,/]  !l 

D5+l,:  =  Ds,i  h  +  [  "3  C5iI  +  "4  Gsl  +  U5  DS!  ]  I] 


(2.52) 


where 


Ul  =  -  (Esj)    '  I  Tl    T2  ] 


-£?,' 

"2  = 

4 

"3  = 

-   C   Tj     T2 

]  ^, 

(*-£./     T4 

-Ed 

t3) 

"4  = 

[  <£M  ^ 

-  r 

Es 

1 

"5  = 

U4  2/2 

(2.53) 


and  the  (s+  l,t)  prediction  error  powers  are  recursively  calulated  as  follows 


r-,  +  u-,  r, 


2  c4 


rid 

Es+\,t  =  Eit  +  E  Tl    T2  3U3r+  "4^3  +  U5T4 


(2.54) 


The  matrices  I,  and  I2  are  of  dimension  (s  +  t  +  2  x  s  +  /  +  3).  They  are  introduced  to 
provide  symmetry  to  the  matrix  algebra  and  preserve  initial  condition  calculations.  We 
design  the  matrices  to  perform  the  following  operations 
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[  wj  ...  wJ+1  co5+2  ...  us+[+2  3  I,  =  [  ft),  ...  ojs+]  0  wJ+2  ...  coJ+f+2  ]  (2.55) 

[  oj)  ...  <wJ+1  w5+2  ...  w5+r+2  :  I2  =  [  0  w,  ...  coJ+1  0  wJ+2  ...  cos+[+]  ]  (2.56) 

The  values  t,  through  r4  can  be  obtained  through  the  correlation  data  and  forward  and 
backward  estimation  parameters.  We  express  them  mathematically  as 

Cx,  T2:  =  £[^r(^-l)vv(A)] 

t3  =  -£[^-5-1)^(A-)]  (2.57) 

T4  =  EZy(k-S-l)dSit(k)l 

The  MA-type  recursive  formula  for  the  forward  estimation  parameter  matrix  and  back- 
ward estimation  parameter  vectors  is  obtained  in  a  similar  manner.  The  recursive  for- 
mula is  given  by 

Vr,r+]  =  ^s,t  ^3  +  nl  Q*,/ 14 

Dv+1  =  [DJjr  +  «2GJ>r]l3  (2.5S) 

Gs,!+\  =  Gs,  h+l  n3  CSJ  +  «4  Ds[  +  A75  G5  f  ]  I3 


where 


"i  =  -  (Eit)']  [  f,  r'2  ] 


/72 

= 

n3 

= 

"[-'] 

T'2 

]  Ev 

(£?/ 

*4" 

-£?,t' 

3) 

(2.59) 


[  (£f,,")2  -  £j,  <,  ] 
«5  =  nA  n2 

and  the  (s,t+  1)  prediction  error  powers  are  calculated  using  the  following  recursive  for- 
mulas 

Es,,+i=Ev+nr[>'i  ^'ii 

^s,t+\  =  T  3  +  n2  T '4 

£?,r-H  =  £V  +  [  T'l   T'2  ]  n3  +  "4  T'3  +  «5  t'4 
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where  I3  and  I4  are  (s  + 1  +  2  x  s  +  t  +  3)  dimensional  matrices  which  we  have  designed 
to  perform  the  following  operations 

[  ctfj  ...  a>5+1   cos+2  ...  cos+l+2  ]  I3  =  [  <y,  ...  w5+1   u)s+2  ...  cos+[+2  0  ]  (2.61) 

[  co,  ...  o>5+1   ws+2  ...  (os+[+2  ]  I4  =  [  0  «!  ...  os  0  co5+2 ...  a)s+!+2  ]  (2.62) 

The  values  of  x',  through  t'4  are  calculated  using  correlation  data  in  conjunction  with 
current  forward  and  backward  parameter  estimation  values.  We  express  these  quantities 
mathematically  as 

[t'j  t'2]  =  -£[^;(A-1)vJi;(A)] 

t'3  =  -E  [  x(k  -  t  -  1)  dsJk)  ]  (2.63) 

T'4=Elx[k-t-\)sJk)2 

It  is  interesting  to  note  that  the  MA-type  recursive  formula  is  the  complimentary  form 
of  the  AR-type  formula  and  that  the  two  are  identical  if  the  variables  associated  with  the 
input  signal  x(k)  and  the  variables  associated  with  the  ouput  signal y(k)  are  interchanged. 
That  is.  we  replace  v(A:).  GJf  and  gSJ(k)  with  —  x(k),  D5t  and  d5,(k)  and  vice  versa. 
1.     Experimental  Results 

The  ARMA  parameter  estimation  algorithm  of  [Ref.  4:  pp.  619-621]  based  on 
the  recursive  formulas  of  equations  (2.52)  and  (2.58)  was  implemented  using  the  Fortran 
program  found  in  Appendix  A.  This  program  calls  subroutines  which  compute  the 
ARMA  model  parameters  as  the  AR  order  is  increased  by  one  and  as  the  MA  order  is 
increased  by  one.  These  subroutines  are  shown  in  Appendix  B  and  Appendix  C.  respec- 
tively. In  the  main  program,  an  input  data  sequence  of  white  Gaussian  noise  is  passed 
through  a  known  reference  model  producing  an  ouput  data  sequence.  We  obtain 
autocorrelation  and  crosscorrelation  data  from  these  input  and  output  sequences.  The 
correlation  data  is  used  to  calculate  initial  values  of  the  error  powers  for  x  and  y  as  well 
as  r,  through  t4  and  t\  through  r'4  .  Next  we  obtain  estimates  of  the  reference  model  by 
employing  the  recursive  formulas  (2.4S)  through  (2.50),  (2.52)  and  (2.58).  Several  refer- 
ence models  were  estimated  beginning  with  a  strictly  AR  process  of  order  5  =  4  having 
as  its  transfer  function 

H(z) ~ ~-J ~ it  (2.64) 

l-0.2z   '+0.62z      -  0.152  z  3  +  0.3016  z 
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The  actual  values  of  the  reference  model  parameters  and  the  ARMA  model  parameters 
which  estimate  this  reference  model  are  listed  below 


ACTUAL  ESTIMATED 

AR:     0.2000  0.2005831 

-0.6200  -0.6207655 

0.1520  0.1527565 

-0.3016  -0.3020376 

MA:     1.0000  1.0001506 

We  next  consider  a  second  reference  model  with  MA  order  /  =  2  and  AR  order  5  =  3 
having  transfer  function 

0.5-040,-  +  0.89,-*  (165) 

1  -  0.20  z     -  0.25  z     +  0.05  z 

The  true  reference  model  parameters  and  ARMA  model  parameter  estimates  are  shown 
to  be 

ACTUAL  ESTIMATED 

AR:  0.2000  0.1993060 

0.2500  0.2496567 

-0.0500  -0.0491961 

MA:  0.5000  0.5002602 

-0.4000  -0.3997071 

0.8900  0.8S94749 

A  third  example  with  MA  order  /  =  2  and  AR  order  s  =  4  having  transfer  function 

1  +  0  2  --1  —  0  99  -~2 

H(z)  = 1  +  "V2 ^^ (2.66) 

l-0.2z   '+0.62z     -  0.152  z  3  + 0.3016  z  4 
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was  considered  for  which  we  obtained  the  following  actual  and  estimated  reference 
model  parameter  values 


AR: 


MA: 


ACTUAL 

ESTIMATED 

0.2000 

0.2011805 

-0.6200 

-0.6223S03 

0.1520 

0.1534197 

-0.3016 

-0.3036S23 

1.0000 

0.999763S 

0.2000 

0.1998342 

-0.9900 

-0.9SS6S52 

We  consider  as  a  final  example  the  reference  model  of  AR  order  and  MA  order  5  =  3  and 
t  —  3,  respectively,  with  specific  transfer  function 


H(z) 


0.5  _  q.95  r~'  +  1.33  z~2  -  Q.979  r~3 
1  +  1.69  -_1  -  0.962  z~2  +  0.2  z"3 


(2.67) 


The  actual  and  estimated  ARM  A  parameters  are 


ACTUAL 

ESTIMATED 

AR:    -1.6900 

-1.6981325 

0.9620 

0.969099S 

-0.2000 

-0.2018440 

MA:     0.5000 

0.4995653 

-0.9500 

-0.9553509 

1.3300 

1.3346767 

-0.9790 

-0.9S64S98 

The  above  examples  demonstrate  the  validity  of  the  parameter  estimation  algo- 
rithm of  [Ref.  4:  pp.  619-621].  Many  reference  models  were  estimated  using  this  algo- 
rithm, including  pure  MA  processes,  for  which  accurate  estimates  were  obtained. 
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D.     LATTICE  STRUCTURE 

In  section  C  we  developed  expressions  for  the  forward  and  backward  prediction  er- 
rors, namely,  those  of  equations  (2.31),  (2.32)  and  (2.33).  From  these  prediction  error 
equations  we  can  design  elementary  AR,  MA  and  AR.V1A  lattice  structures  or  sections. 
Each  elementary  section  satisfies  the  orthogonal  conditions  as  listed  in  equation  (2.37). 
From  the  prediction  error  recursive  formulae,  equations  (2.31),  (2.32)  and  (2.33),  we 
construct  the  AR-type  elementary  lattice  section  as  follows.  Consider  the  following  data 
set  of  order  (5+  1,0  consisting  of  input  and  output  data  elements, 

hs+u(k)  l{=  ly(k)  ...y(k  -  s)     -x(k)  ...  -*(*  -  /  +  1)  -x{k  -  t)  ]  (2.68) 

WA')  l2  =  CK*  "  1)  •••  Ak  ~  s  ~  1)     "*(*  -  0  .«  -*(*  "')<>]  (2.69) 

where  If  and  IJ  are  the  transposes  of  the  matrices  Ij  and  \2  defined  in  equations  (2.55) 
and  (2.56).  We  obtain  a  recursive  relationship  between  the  forward  prediction  errors 
\(k)  of  order  (s  +  \.t)  and  order  (s,t)  by  substituting  equations  (2.52),  (2.68)  and  (2.69) 
in  equation  (2.31)  such  that 

=  hs+]Jk)[\lcl[  +  l2Dl!u]]  (2.70) 

=  \t(k)  +  dStl{k-  l)u, 

The  backward  prediction  error  recursions  are  obtained  in  a  similar  manner  and  the 
AR-type  error  recursions  are 

vs+u(k)  =  vs,t(k)  ~  u\  ds>lk  -  1) 

4i#W,,(*)  +  ^4,f(*-i) 
Ss+\,t(k)  =  ssJk)  -  ui  dsAk) 

w*) = <w*  - 1) + [  «f  ^  ]  vv(A-)r-  u4g5+lit(k) 

where  u,  =  E  u\  u\  ~J  and  u3  =  [  u\  u\  ]■  The  AR-type  elementary  lattice  inverse  section 
based  on  these  error  recursions  is  shown  in  Fieure  5. 
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v     f(k) 

s,t 


g   t(k) 

s,t 


vY   (k) 
s,t 


ds,t^ 


v       .     (k) 
■►     s+1,t 


g        (k) 

-►    s+1,t 


v  (k) 

**     s+1,t 


d      1  <(k> 
-►     s+1,l 


Figure  5.      AR-type  elementary  lattice  inverse  section 

We  design  the  MA-type  elementary  lattice  inverse  section  in  a  similar  manner  using  the 
following  representation  of  the  data  set  of  order  (s,t  +  1) 


'j,;+] 


k)  l[=  [ y(k)y{k  -  1)  ...y(k  -  s)     -x(k)  ...  -.x(k  -  t  +  1)    -jc(A  -  i)  ] 


h, ..,  (AM  Ij  =  [  jl  A  -  1)  ...y(A  -  5)  0    -x{k  -  1)  ...  -.v(A  -  ;  -  1 )  ] 


(2.72) 


and  by  substituting  equations  (2. 58)  and  (2.72)  into  equation  (2.31).  we  obtain  the  for- 
ward prediction  error  recursion  as  the  MA  order  is  increased  by  one,  namely, 


(2.73) 


The  backward  prediction  error  relationships  are  obtained  in  a  similar  manner  and  are 
iiiven  bv 


4.r+l(*)"<W*)-w2&.i(*) 

S«+l(A)  =  &,/*  ~  1)  -  [  "3*   »3  ]  vsAk)T-  /74  ^.H-lW 


(2.74) 


The  MA-type  elementary  lattice  inverse  section  based  on  these  error  recursions  is  shown 
in  Ficure  6. 


20 


v       (k) 


g      (k) 
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Vs,t(k) 


d       Ik) 
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>Vs,t+lW 


^9     t    A*) 
^    s,t+1 


s,t+1 


(k) 


d  (k) 

*■    s,t+1 


Figure  6.      MA-type  elementary  inverse  lattice  section 


We  now  construct  the  ARMA  elementary  lattice  section  from  the  AR  and  MA 
prediction  error  recursions.  Assuming  that  the  prediction  errors  are  known  for  a  given 
model  order  (s,t),  the  (s+  l.t)  prediction  errors  can  be  calculated.  These  prediction  errors 
of  order  (s+  l.t)  are  then  updated  as  the  MA  order  increases  by  one  resulting  in  a  pre- 
diction error  of  order  (s-r  l,t+  1).  We  consider  the  forward  prediction  error  for  x  as  the 
AR  order  is  increased  by  one.  specifically 


&,,(*)  =  i£(A)-i*f<arw  (A-  1) 


2.75) 


Now.  v;.! ,  (A)  becomes  the  current  value  of  the  forward  prediction  error  for  .y  and  when 
we  calculate  the  (s+  l,t+  I)  forward  prediction  error  we  have  from  equation  (2.73) 


v*+u+1  (A)  =  vf+u  (A)  +  n?gs+u  (A  -  1) 


(2.76) 


Equation  (2.76)  can  be  expressed  in  terms  of  the  (s.t)  forward  prediction  errors  of  x  by 
making  appropriate  substitutions  for  v;_lf(A)  and  g,,u,(k-  1).  That  is,  we  substitute 
v'.u,(k)  and  gs.lt,(k-  1)  of  equation  (2.71)  in  equation  (2.76)  to  obtain  the  (s+  l,t+  1) 
forward  prediction  errors  for  .x,  namely, 

v?+u+I  (A)  =  i£  (A)  -  «f  dStt  (A  -  1)  +  nf [  gs<l  (A  -  1)  -  u2  dSJL  (A  -  1)  ]         (2.77) 
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Grouping  the  terms  we  obtain 

W,  (k)  =  <,  (A)  -  («f  +  «f  «j)  </w  (A-  -  1)  +  n*gs,t  (*  "  1)  (2.78) 

The  forward  prediction  error  recursion  for  y  of  order  (s+  l.t+  1)  is  obtained  in  a  similar 
manner.  We  begin  with  the  (s+  l,t)  order  update  of  the  prediction  error  and  after  it  is 
computed,  update  the  MA  order.    Specifically,  we  have 

vy5+u(k)  =  ^,(k)  +  uvlds,(k-l)  (2.79) 

and  from  equation  (2.73) 

v?+1>r+]  (k)  =  4U  (A-)  -  nylgs+u  (*  "  0  (2-8°) 

Substituting  (2.71)  for  v>_lr(A)  and  g^ltt(k  —  1)  in  equation  (2.80)  then  grouping  terms 
we  obtain  the  (s+l,t+I)  forward  prediction  error  recursion  for y, 

&i,,+i  (*)  =  <r  (A)  +  («f  +  «f  «2)  4>f  (A  -  1)  -  «f  &tl  (A  -  1)  (2.S1) 

The  (s+  l,t+  1)  backward  prediction  errors  for  x  and  y  are  derived  in  a  similar  manner 
and  are  given  by, 

ft+i.r+1  (A')  =  «m  (*  -  1)  +  ("3T  +  »4  "?)  "J  (*)  -  («3  +  >h  4)  vji  (k)  ^ 

4+i, ;+i  (*)  =  ds>l  (A  -  1)  -  w3x  v*,  (A)  4- 1?3  vj,  (A) 

The  ARM  A  elementary  lattice  inverse  section  is  shown  in  Figure  7  where  the  coefficients 
are  related  to  the  prediction  error  recursions  by  the  following 

1       /   x   ,      x      \  1  x  1        /    v    .      y      \  1  v  /  ->  c  -  \ 

vv-j  =  (ux  +  /7,  w2),     iv2  =  nx ,     vv3  =  (tq  +  n\  u2),     w4  =  n\  (2.5j) 

w\  =  (n*  +  w4  uf),     w\  =  (a?3  +  «4  «£),     ivj  =  w*,     vvg1  =  wf  (2.84) 
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Figure  7.      ARMA  elementary  lattice  section 

We  see  from  Figure  7  that  each  elementary  ARMA  lattice  inverse  section  contains  eight 
coefficients. 

From  the  AR,  MA  and  ARMA  elementan."  lattice  inverse  sections,  we  can  obtain 
synthesis  lattice  structures.  These  structures  provide  a  means  of  working  with  lattice  re- 
alizations as  linear  filters.  The  resulting  AR,  MA  and  ARMA  elementan'  synthesis 
lattice  filters  are  shown  in  Figures  S  and  9  respectively. 
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Figure  8..     Top:  AR  elementary  lattice  filter.  Bottom:  MA  elementary  lattice  section. 
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Figure  9.      ARMA  elementary  lattice  filter 

Summarizing,  in  this  chapter  we  have  reviewed  the  Mullis- Roberts  criterion,  intro- 
duced the  ARMA  parameter  estimation  as  a  generalized  Mullis-Roberts  criterion  and 
obtained  analysis  and  synthesis  forms  of  lattice  structures.  We  notice  that  each  ARMA 
elementary  lattice  section  consists  o[  eight  reflection  parameters  and  the  calculation  of 
these  parameters  requires  the  autocorrelation  and  crosscorrelation  information  as  ob- 
tained from  the  input  output  data  of  the  reference  model.  Also,  we  obtained  a  set  of 
equations  relating  the  final  model  estimation  parameters  and  the  prediction  error  model 
parameters  which  inturn  are  obtained  from  the  eight  lattice  parameters. 
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III.     ADAPTIVE  LATTICE  ALGORITHM 

A.     LEAST  MEAN  SQUARE  ALGORITHM 

The  study  and  design  of  adaptive  filters  is  known  to  be  a  very  important  part  of 
statistical  signal  processing.  Many  adaptive  algorithms  have  been  developed  to  support 
the  application  of  adaptive  filtering  in  communications  and  control  [Ref.  12].  An  adap- 
tive filter  is  characterized  by  the  ability  of  its  filter  coefficients  to  adjust  (self-optimize) 
automatically  and  yield  an  optimum  filter  design.  Two  processes  occur  within  an  adap- 
tive filter,  namely,  the  adaptation  and  the  filtering  processes.  During  the  filtering  process 
a  desired  signal  is  applied  to  an  adaptive  algorithm  as  a  reference  for  adjusting  the  filter 
coefficients.  Figure  10  shows  a  block  diagram  of  the  adaptive  modeling  process.  Refer- 
ring to  Figure  10.  let  y(k)  be  the  output  of  the  filter  at  time  k.  By  comparing  the  output 
with  the  desired  signal  d{k)  .  an  error  signal  e(k)  is  generated.  The  adaptive  algorithm 
of  the  filter  uses  this  error  signal  to  generate  corrections  which  are  applied  to  the  filter 
coefficients  such  that  an  optimum  solution  is  obtained.  An  optimization  technique  called 
the  method  of  steepest  descent  provides  an  approach  to  solving  this  problem.  The  pro- 
cedure is  as  follows: 

1.  Assign  initial  values  to  all  filter  coefficients. 

2.  Using  these  initial  values,  compute  the  gradient  vector,  whose  individual  elements 
equal  the  first  derivatives  of  the  mean-squared  error  with  respect  to  the  filter  coef- 
ficients. 

3.  Compute  values  for  the  filter  coefficients  by  changing  the  initial  values  in  the  di- 
rection opposite  that  of  the  gradient. 

4.  Return  to  step  2  and  repeat  the  procedure. 

There  is,  however,  a  limitation  to  this  procedure.  The  steepest  descent  algorithm  requires 
exact  measurements  of  the  gradient  vector  at  each  iteration  which,  in  practice,  is  not 
possible.  Therefore,  the  gradient  vector  must  be  estimated  and  consequently,  errors  are 
introduced.  An  algorithm  is  required  which  computes  the  gradient  from  the  available 
data.  The  least  mean  square  (LMS)  algorithm,  developed  by  Widrow  and  Hoff,  is  widely 
used  and  is  very  convenient  to  implement  in  real  time  hardware  [Ref.  13:  pp.  96-104]. 
Let  y(k)  be  the  output  of  the  filter  and  d{k)  the  desired  signal  at  time  k  as  shown  in 
Figure  10.  We  compute  the  error  by  taking  the  difference  between  these  two  signals, 
namelv. 
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Figure   10.      Adaptive  modeling  block  diagram 


c(k)  =  d{k)  -  y{k) 


(3.1) 


The  value  of  the  mean-squared  error  is  the  expected  value  of  the  error  squared, 
•£[  e2{k)  ]  an^  the  gradient  vector,  V  (A)  ,  is  the  first  derivative  of  the  mean-squared  er- 
ror. The  eradient  vector  is  siven  bv 


^=i£[^]=2^^lir^ 


(*) 


(3.2) 


"where  w(k)  is  the  time  dependent  filter  coefficient  vector.  The  recursion  for  the  filter 
coefficient  which  changes  the  old  value  in  the  direction  opposite  to  that  of  the  gradient 
is  then  given  by, 


yy(k)  =  ys(k-l)  +  j-  vl  -V  (A)  ] 


=  yy(k-\)-ne(k) 


6\\{k) 


e(k) 


(3.3) 


where  \\(k)  is  the  filter  coefficient  vector  estimate  at  the  k'h  iteration,  w(k  —  1)  is  the  past 
filter  coefficient  vector  estimate,  ^  is  the  convergence  (gain)  constant,  e(k)  is  the  error 
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sienal  at  the  k"'  iteration  and    .  '       e(k)  is  the  instantaneous  gradient.    The  implemen- 
c  c\\{k)  fc  h 

tation  of  this  algorithm  proceeds  as  follows: 

1.  Assign  initial  values  to  the  filter  coefficients. 

2.  Compute  the  value  for  the  error  signal  e(k). 

3.  Calculate  the  updated  estimate  of  the  filter  coefficients  using  the  instantaneous 
gradient. 

4.  Increment  the  time  index  by  one  and  return  to  step  1. 

Convergence  properties  of  the  LMS  algorithm  are  well  documented  within  the  literature. 
The  choice  of  a  gain  constant  n  is  arbitrary  however,  theoretical  bounds  have  been  de- 
rived for  n  ,  given  by  [Ref.  14:  pp.  101-106], 

0  <n  <~2—^  Jrrl    1  (3.4) 

'max  lT\-Kxx-i 

where  /max  is  the  maximum  eigenvalue  of  the  input  autocorrelation  matrix,  R„.  and 
where  Tr  [  R„  ]  is  the  trace  of  the  matrix  R„. 

B.     DERIVATION  OF  THE  ADAPTIVE  LATTICE  ALGORITHM 

The  adaptive  lattice  algorithm  developed  in  this  thesis  uses  concepts  of  the  LMS 
algorithm  discussed  in  section  A  and  applies  them  to  the  ARMA  digital  lattice  filter 
proposed  in  Chapter  II.  Consider  the  ARMA  digital  lattice  filter  of  Figure  11,  which 
consists  of  two  cascaded  elementary  lattice  sections.  The  filter  coefficients  (weights)  are 
defined  such  that  wf  represents  the  i"'  lattice  coefficient  at  stage  m  of  the  lattice  structure. 
In  this  figure  we  have  a  two  stage  lattice  and  there  are  eight  coefficients  per  elementary 
lattice  section.  The  output,  y{k),  of  the  lattice  filter  can  be  determined  from 

y(k)  =  ftk)  +  w\  e^k  -  1)  -  w\  4o(k  -  1)  (3.5) 

Forward  errors  at  a  given  stage  m  of  the  lattice  filter  are  defined  as, 


±  (*)  -  <L,  w  +  ""+'  <(*-»)-  <+'  <  (*  -  » 


(3.6) 


and  the  backward  errors  for  any  given  stage  m  are, 


<  w  -  <-,  (*  - «> + <  eL  w  -  «f  '£-  w 

<  (*)  =  <C,  (*-')-  <  «/„-,  (*)  +  «f  «L  <*) 


(3.7) 
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Two  stage  ARMA  lattice  digital  filter. 
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where,  in  Figure   11.  m=  1.2  and  e}  =  x(k)  and  e}  =y(k).  The  terminal  condition  is 
ef  (k)  =  b0  e)  (k);    m  —  1  .  To  begin  with,  let  l\  equal  unity.  The  initial  conditions  are 
e\  (k  —  1)  =  0  and  el  (k  —  \)  =  0.  As  with  the  LMS  algorithm,  we  form  an  error  between 
a  desired  signal  d(k)  and  the  output  signal  y(k)  such  that, 

e(k)  =  d{k)-$(k)  (3.8) 

The  instantaneous  gradient  according  to  eq  (3.2)  is  then, 

V(A)  =  2  e(k)  -^y  [  d{k)  -y(k)  ]  (3.9) 

Since  the  desired  signal,  d(k),  is  not  a  function  of  the  filter  coefficients,  equation  (3.9) 
reduces  to. 

V(k)  =  2e(k)-^—  [  -#*)]  (3.10) 

where  the  quantity    „  "  ,      f    —  y(A)  1  is  referred  to  as  the  gradient  estimator.  This  sra- 

1  '    cv\(  k)    L       *        J 

dient  estimator  must  be  computed  for  each  filter  coefficient  within  the  lattice  structure. 
The  filter  coefficients  are  then  updated  using  the  respective  gradient  estimators.  That  is. 

we  need  to  compute. 

V(A)  =  -f—  y  (A)  for/=  1.2, ...  .8  and;  =1,2 M  (3.11) 

3  wj 

where  M  is  the  number  of  stages  in  the  ARVIA  lattice  filter.  From  equation  (3.5).  the 
gradient  estimates  are  given  bv 


rv(A)           d${k)       dwl      _                     ,    dexh{k-\)       dw] 
-^-r-    =  '    .      + 7-  ex  (k  -  1)  +  wl : L  4  (k-\) 

GUT;  C\Vj  OWj  CWj  CW; 

\  y(l       n    (3-12) 
,    c^(A-l) 


—  vv- 


Clt; 


Let  i//(A)  represent  the  partial  derivatives  of  the  output  y(k)  with  respect  to  the  filter  co- 
efficients and  4>{k)  represent  the  partial  derivatives  of  the  errors  with  respect  to  the  filter 
coefficients.  Using  this  representation  we  can  re-write  equation  (3.12)  as  follows, 
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hj  (*)  =  </>/,  ij  (A)  +  6\\  <  (A  -  1)  +  wl  Vb,ij  (k  -  1)  -  Sj{  <  (A  -  1) 

i     v  (3.13) 

where  <5J{,  (5j;,  are  kronecker  deltas  whose  value  is  one  if  and  only  if 
/  =  4  and  j  =  1  or  /'  =  3  and  y  =  1  respectively.  We  compute  \ptJ  (k)  ,  by  obtaining  re- 
cursive relations  which  calculate  the  partial  derivative  of  the  forward  and  backward  er- 
rors with  respect  to  the  filter  coefficients.  These  relations  are  obtained  from  equations 
(3.6)  and  (3.7)  by  taking  the  partial  derivatives  with  respect  to  the  filter  coefficients, 
namely, 

$. « (*) = c  u  w + «?/  o*-  d + «?  «v,  (/*(*-«)-  yy  <_,  t*  - ») 
€  u  (*) = c « w + «sr,u  <  (*-d+ *™+'  *jl  v  (*  - 1)  -  «+l,;  4,  (k  - 1) 

-w3m+1  4>yb  u{k-  1) 


<. y  W  -  «L, u  (*  -  D  +  *s  / £_,  (*)  +  <  <%_, u (*) "  *"  / 4L.  (*) 


■w?  #_,/;(*) 


(3.14) 


These  recursive  relations  possess  a  lattice  structure  similar  to  that  of  Figure  11.  with 
delta  components  injected  at  the  summation  nodes.  They  may  also  be  simplified  by  ex- 
amining individual  terms.  Consider  the  general  equation  for  the  forward  error  in  x,  re- 
peated here  for  continuity, 

<£  w -  eL  w  +  < C (*-«-  *T <-, (* -  o  (3-15) 

The  partial  derivative  with  respect  to  each  filter  coefficient  is  expressed  as, 


.in 

x 


c\vt 


cvv,         „  _, 


7«^x(A-l)-w, 


<-,  (*  - 

1) 

Owl 

1) 

(3-16) 


CVV-  OW[ 
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Since  the  partial  derivative  is  taken  with  respect  to  the  current  filter  coefficient  w;  at  time 
k,  the  partial  derivatives  involving  delay  terms  i.e.,  (A  —  1).  are  set  to  zero.  This  result 
follows  from  the  realistic  assumption  that  eim  (k  —  1)  is  a  function  of  w;  (A  —  1)  but  not 
of  w>t  (k)  .  Also  note  that  wi  (k)  is  a  function  of  w{  (k  —  1)  but  not  vice  versa.  With  these 
simplifications  we  reduce  the  equations  of  (3.14)  to. 


*/.v(*)-*JLv  (*)  +  «?/ <_,<*-  ">-*"/ <_,<*-  1) 

*«(*)-<iL«(*)+*t,w<(*-1)-«f,M<(*-») 

<  « (*)  -  «#  C  (*) +  <  C,.j(*)-*i,w-  <"  C  « (*) 

<  v  (*)  =  -  *tf  <£.,  w  -  «T  C  u  (*)  +  *  <&.,  w  +  *™  C  <j  <*) 


(3.17) 


and  the  gradient  estimator  is. 

yij(k)  =  tfiU(k)  +  d\Jie!o(k-l)-dlJiett)(k-\)  (3.18) 

Although  these  are  valid  recursive  relations,  they  are  difficult  to  implement  in  a  lattice 
algorithm.  The  ultimate  goal  is  the  requirement  to  easily  compute  ^i:{k)  from  the 
available  data.  From  eq  (3.18)  it  is  evident  that  t//,7  (A-)  depends  on  <p}{,}  (A)  which  in  turn 
requires  knowledge  of  4>»f  tJ  (A)  and  <f>*f  u  (A)  but  not  of  <j>i  u  (A)  or  4>*b  fJ  (A).  Therefore  the 
three  equations  necessary  to  compute  the  gradient  estimator  are, 

hj  (*)  ~4%u  (*)  +  sit  <  (^  - 1)  -  *J  /  <  (*  - !) 

€/y(A)  =  C1/yW  +  ^7<_I(^-i)-^/<_1(^-i)  (3.19) 

These    equations    are    dependent    on    the    filter    coefiicients    viv.     /=  1.2.3.4    and 

y=  1.2, M  ,  thereby  reducing  by  one-half  the  number  of  computations  required  for 

<£;  7  (A)  and  <j>f  ,7(A).  A  recursive  relation  is  desired  for  \jji;  (A)  which  does  not  involve 
delta  functions.  Consider  the  four  stage  lattice  filter  with  terminal  condition  b0  equal  to 
unity  such  that  <^V,,-(A)  =  <f>}MlJ -(A)  .  The  procedure  for  computing  \jjtJ  (A)  using  the 
equations  of  (3.19)  is  as  follows: 

1.  Calculate  </>;  ,7(A)  with  m  equal  to  one  and  letting  i  and  j  range  from  one  to  four. 
Repeat  for  m  =  2,3,4. 

2.  Using  the  terminal  condition  and  expression  for  ty  tJ  (A),  calculate  (j>iu{k). 

This  requires  solving  88  equations,  however,  the  result  is  a  very  simple  recursive  formula 
for  if/jj  (A)  .  namely, 
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<{,ij{k)  =  -eybj_l(k-l)         f-1,3 

^ij{k)  =  exbjJk-\)         i-2f4  (J,20) 

The  lattice  coefficients  are  calculated  by  substituting  the  recursive  formula  of  eq  (3.20) 
for  the  gradient  estimator  in  equation  (3.3).  The  adaptive  coefficient  update  equation  is, 

n-j(k)  =  nJi(k-l)-ne(k)ii,iJ(k)  (3.21) 

Since  the  gradient  estimator,  and  therefore  the  gradient,  is  the  same  for  w>  ,  /=  1,3  and 
w;.  i  =  2,4  ,  it  follows  that  w{  =  w{  and  \v2  =  w4.  The  number  o[  filter  coefficients  re- 
quired to  update  the  lattice  filter  is  reduced  to  two,  i.e.,  w[  and  w2.  Furthermore,  from 
the  symmetry  of  the  lattice  structure,  the  following  equalities  between  filter  coefficients 
are  assumed, 


4 

= 

4 

4 

= 

4 

4 

= 

4 

(3.22) 


\vi    =    w. 


Incorporating  these  equalities  with  those  derived  by  the  gradient  estimator  produces  the 
elementary  ARMA  lattice  section  of  Figure  12  where. 

*2  =  4  -  4  =  4  -  rj 

J       J       J       J      I  {'-2^ 

To  prove  that  these  coefficient  reductions  are  valid,  a  computer  generated  solution  using 
the  Fortran  program  of  Appendix  D  was  compared  to  hand  analysis  of  a  second  order 
transfer  function  and  lattice  filter.  The  output  of  the  ARVIA  digital  lattice  filter  was  first 
put  into  difference  equation  form  and  then  compared  to  the  known  transfer  function. 
From  this  comparison,  lattice  coefficients  were  computed.  Details  of  this  analysis  are  as 
follows.  Consider  a  two  stage  ARMA  digital  lattice  filter  comprised  of  the  reduced  ele- 
mentary section  shown  in  Figure  13  and  a  transfer  function  of  the  form. 

bn  +  b,  z~  +  b->  z~ 

l+fl,2         +^2 
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Figure   12.      Simplified  elementary  ARMA  lattice  section. 

The  output  of  the  lattice  filter  can  be  written  in  difference  equation  form  by  earning  out 
the  following  steps:  (i)  start  with  the  output  of  the  lattice  filter  equation  (3.5),  (ii)  sub- 
stitute expressions  for  the  forward  and  backward  errors,  equations  (3.6)  and  (3.7).  re- 
spectively, into  equation  (3.5)  and  (iii)  carry  out  the  algebra.  A  detailed  derivation  of  this 
difference  equation  is  given  in  Appendix  E.  The  difference  equation  in  its  final  form  is 
given  by, 


y{k)  =  x{k)  +  2{  w\  +  w]  w]  +  w\  w\  )x(k  -  1)  +  2  vv;2  x(k  -  2) 
— 2(  v,-,   +  w2  vv2  +  vv-j  u,  )y-(k  —  1)  —  2  wj".v(A'  —  2) 


(3.25) 


which  can  be  written  in  the  transfer  function  form  as, 
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H{z)  = 


1  +  2(  vv2    +   VV,  ;v,    4-   vv2  w2  J  z      +2  vv2  z 


1  +  2(  VV,    +    vv2  n;2    +    ^'i  W]  )  2       +2  vvf 


-2 


(3.26) 


Comparing  the  lattice  filter  transfer  function,  //,(z)  with  the  known  filter  transfer  func- 
tion Hj(z),  produces  the  following  relationships  between  filter  coefficients, 


bx  =  2(  vv2  +  W]  W\   +  w2  w2  ) 


_  o 


1       2 


2(  vvj    +  vv2  w2   +  Wj  vv,  ) 


b2  =  2  w~ 


(3.27) 


a2 


vv. 


Solving  for  vv?  and  vvf  in  terms  of  the  known  transfer  function  coefficients  b2  and  a2  and 
then  substituting  these  results  into  the  expressions  for  bx  and  au  respectively,  yields  the 
following. 


bx  =  a2  vv,    +  (  2  +  b2 )  vv2 

a,  =  (  2  +  a2 )  vv,    +  b2  w2 


(3.28) 


or  in  matrix  form, 


2  +  lh 


2  + a-, 


l 
vv, 

1 
w0 


(3.29) 


and  solving  for  the  lattice  coefficients,  we  have 


v.', 


vv. 


1 


a2b2  -  {2  +  b2){2  +  a2) 


b2  -  (2  +  b 

—  (2  +  a2)  a2 


(3.30) 


and 


Wt   = 


VV,    = 


(3.31) 


Now  that  a  method  of  converting  between  lattice  and  transfer  function  coefficients  for 
a  second  order  system  has  been  established,  we  consider  the  specific  transfer  function 
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1   -  0.8  z~l  +   1.73  r'2  . 

7  1  -  0.89  2   '  +  0.25  z  2 

where 

/?0=1.0      £,  =  -0.80     62  =  1.78 
al  =  -0.89      a2  =  0.25 

From  equations  (3.30)  and  (3.31)  the  lattice  coefficients  are  calculated  as. 

w\  =  0.125,       w\  =  0.890,       wj  =    -0.240719,       wj  =    -0.195719 

Values  for  the  steady-state  lattice  coefficients  were  computed  using  the  Fortran  program 
in  Appendix  D  and  are  shown  below.  Convergence  aspects  of  both  the  lattice  coeffi- 
cients and  output  error  are  shown  in  Figure  13. 

w\  =0.124982,     w\  =0.890003,     w\  =-0.240710,     w\  =-0.195711 

these  results  confirm  the  validity  of  the  derived  adaptive  lattice  algorithm  and  the  design 
of  a  new  elementary  lattice  section  shown  in  Figure  12. 

The  current  adaptive  lattice  algorithm  assumes  that  the  terminal  condition  is  unity. 
This  is  generally  not  the  case  in  practice.  We  now  extend  this  adaptive  algorithm  to  the 
more  general  case  where  the  terminal  condition  is  an  arbitrary  constant.  The  recursive 
relation  which  updates  bD  is  similar  to  those  which  update  the  other  lattice  filter  coeffi- 
cients. The  update  equation  for  b0  is  given  by 

ce(k) 
bQ(k)  =  b0(k  -  1)  -  n  e(k)  -j-—  (3.33) 

de{k) 

The  sradient  estimator  — — — ,  is  calculated  usins  equations  (3.5),  (3.8)  and  the  fact  that 

cb0{k) 

the  desired  signal  dfk)  is  not  dependent  on  b0.  The  gradient  estimator  for  b0  is  written 
as, 


l 


8y(k)       SeJ^k)    |    dw\      X/i       ^         ,    6exb(t(k  -  I)        dw\ 

) 


+  ~, —  ^blb  -  1)  +  Wa  n— -tt^  eUk  -  1) 

cb0  6b0  dbQ      ^  4  cbQ  cb0      °°K 


(3.34) 


—  Wi 


cb0 
Since  the  partial  derivative  is  taken  with  respect  to  b0  at  time  k,  this  reduces  to, 
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cyjk)        d^ 
cb0  dbQ 


(3.35) 


Similarly, 

eL  W  =  *L  (*)  +  vv"+1  <  (*  "  ])  "  M'3W+1  <  (*  -  1)  (3.36) 

and  the  partial  derivative  with  respect  to  b0  is 


560  <?6C 


(3.37) 


The  terminal  condition  is. 


&(*)-*>&(*)  (3.3S) 


Taking  the  partial  derivative  of  equation  (3.38)  with  respect  to  £0  yields  e}  (k),  and  the 
recursive  equation  to  update  ba  (k)  becomes, 

b0(k)  =  b0(k-l)-fie{k)e?M(k)  (3.39) 

The  gradient  estimators  for  the  lattice  filter  coefficients  are  scaled  by  the  arbitrary  con- 
stant b0,  since  at  the  terminal  condition  </>>'  ,  (k)  =  b0  4>xf  u  (A)  and  b0  is  propagated 
through  the  calculations.  The  gradient  estimators  become, 

^lj(k)  =  -b0e>b     {k-l)     i=1.3 

x  (3-40) 

^iJ(k)  =  b04_i(k-\)     i-2,4 

To  test  this  more  general  adaptive  lattice  algorithm  the  output  of  a  known  transfer 
function  with  b0  equal  to  0.5  was  compared  to  the  output  of  the  ARM  A  digital  lattice 
filter.  The  second  order  transfer  function  used  was, 

„,.       0.5  -0.4--'+  0.89  z'2  r,u 

Hj{z)  = —  (3.41) 

J  1  -0.89z   '+0.25  2  z 

The  computer  generated  steady  state  lattice  coefficients  are  given  below  and  convergence 
aspects  shown  in  Figure  14. 

b0  =  0.499946,     w\  =  -0. 120428.     vv22  =  0.593237,     vv,1  =  -0.447423,     w\  =  0. 166320 
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Figure   14.      Top:  Lattice  coefficients.  Bottom:  Output  error. 
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In  order  to  maintain  the  same  adaptive  time  constant  and  misadjustment  at  each 
stage  in  the  lattice,  the  convergence  constant  is  normalized  by  the  power  level  at  each 
stage  [Ref.  15].   Therefore,  we  can  write  equations  (3.21)  and  (3.39)  as, 


J.UA^J."  P 


mi?  (A)  =  ^  (A  -  1)  -  -f—  e(k)  iPiJ  (k) 
a)  (A) 

b0(k)  =  b0(k-l)--£-e(k)el(k) 


(3.42) 


where  p  is  the  convergence  constant  and  o)  (k)  and  y2  (A)  are  estimates  of  the  power  at 
the/*  stage  for  w;  and  b0  respectively  and  computed  as  follows: 

a)  (A)  =  p  o2  (A-  1)  +  (1  -  p)tfj{k) 

y2  (A)  =  p  y2  (A  -  1)  +  (1  -  p)  [  £  (A)  ]2  (3'4jj 

Writing  equation  (3.42)  using  the  notation  adopted  for  the  reduced  elementary  ARMA 
lattice  section  we  obtain 

rj  (A)  =  rj  (A  -  1)  -  -£-  e(k)  exb     (A  -  1) 
oj  (A) 

kj(k)-kj{k-l) j—  e{k)4    (A-l)  f344) 

ffy   (A) 

^0  (A)  -  b0  (A  -  1)  -  ^^  c(A)  «/  (A) 

y  (*) 

In  the  above  equations  p  is  a  weighting  parameter,  0  <  p  <  1,  which  distributes  the 
amount  of  weight  given  the  past  power  level  or  current  sample.  Normalized  convergence 
constants  are  used  in  all  examples  of  this  thesis.  The  adaptive  lattice  algorithm  is  sum- 
marized in  Table  1 . 

In  summary,  we  have  derived  an  adaptive  algorithm  based  on  the  LMS  theory  of 
adaptive  coefficient  computation.  This  new  adaptive  algorithm  easily  updates  the  lattice 
coefficients  by  using  available  data.  The  original  requirement  to  update  eight  coefficients 
of  an  elementary  ARMA  lattice  section  was  reduced  to  updating  only  two  coefficients 
and  still  being  able  to  describe  the  lattice.  The  algorithm  is  general  in  that  it  applies  to 
systems  whose  terminal  condition  is  an  arbitrary  constant.  The  validity  of  this  algorithm 
was  demonstrated  through  comparisons  between  hand  analysis  and  computer  simu- 
lation. In  the  next  chapter,  we  further  demonstrate  the  convergence  of  this  algorithm. 
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Table   1.     SUMMARY  OF  ADAPTIVE  LATTICE  ALGORITHM 

With  given  initial  conditions,  e}Q  =  x(k),  e-}0  =y{k),  e*0  (k  —  1)  =  e?,o  (A  —  1)  =  0  and  all 
lattice  coefficients  zero, 

Step  1:      for  k=  l,m  compute 

«/.«  -  eL  (*)  +  <  <-,  (*  - 1)  -  «T  <-,  <*  -  " 

with  output  ej  (k) 
Step  2:      for  k=  l.m  compute 

<  (*)  =  <-,  (*-d+ »f £.,  (*)  -  <'  <L,  (*) 
«i  (*)  =  <_,  (*-»)-  -"'  <•/„.,  (*)  +  "f  C  (*) 

Step  3:      Update  coefficients 

rj  (A)  =  Tj  (k  -  1)  -  -A-  *(*)  rf  .  (*  -  1) 

a;  (A) 

kj  (A)  =  kj  (A  -  1 )  -  -f—  e(k)  4    (A  -  1) 
a]  (k) 

b0(k)  =  b0(k  -  1)  -  -£-  t-(A)  «/  (A) 

y  (*) 

Step  4:      Repeat  for  next  iteration  i.e.  return  to  step  1. 
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IV.     EXPERIMENTAL  RESULTS 

The  adaptive  lattice  algorithm  derived  in  Chapter  III  is  now  computer  simulated  to 
study  its  convergence  performance.  The  system  identification  mode  of  adaptive  filtering 
is  considered  for  this  purpose.  Figure  15  shows  a  general  system  identification  config- 
uration. The  systems  considered  are  time-invariant  and  linear.  Notice  that  we  apply  the 
same  input,  white  noise  in  general,  to  both  the  reference  system  and  the  adaptive  lattice 
filter  which  is  modeling  the  system.  The  criterion  in  this  configuration  is  to  minimize  the 
mean-squared  error  between  the  system  and  filter  outputs.  Thus,  in  this  context,  the 
adaptive  algorithm  continuously  updates  the  lattice  filter  parameters  in  order  to  mini- 
mize the  mean-squared  error. 

The  adaptive  algorithm  is  realized  as  summarized  in  Table  1.  As  we  mentioned  in 
Chapter  III,  the  two  important  parameters  of  the  algorithm  are  the  adaptation  constant 
M  and  the  weighting  constant  p  .  In  what  follows,  we  shall  consider  convergence  studies 
of  both  second  and  third  order  reference  systems  (fixed  filter  transfer  functions).  Con- 
sider the  following  reference  system  with  transfer  function, 


fl/tf 


1  +  0,2  r"1  -  0.35  z~2 
1  -  1.4  z-1+  0.85  z~2 


This  system  has  complex  poles  and  simple  zeros  located  at  z  =  (0.7  ±y'0.6)  and 
z  =  0.5,  —0.7,  respectively.  Using  a  convergence  constant  n  =  0.01  and  power  level 
weighting  factor  p  =  0.45,  the  adaptive  ARMA  digital  lattice  filter  which  models  the 
above  system  has  the  following  steady-state  lattice  parameters. 

terminal  condition    b0  =  0.999202 

lattice  coefficients    r,1  =  0.3525SO 

r~=  -0.174355 

k\  -  -0.44S228 

Af=  0.425060 

Convergence  properties  of  the  lattice  coefficients  and  error  are  shown  in  Figure  16.  The 
mean-squared  error  was  minimized  after  approximately  1700  iterations  at  which  time  the 
lattice  coefficients  reached  their  steadv-state  values.  When  the  value  of  the  convergence 
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Figure   15.      Block  diagram  of  system  indentification  modeling 


constant  p.  was  modestly  increased,  convergence  degraded  rapidly.  Also,  when  the 
weighting  factor  p  was  increased,  convergence  deteriorated  quickly.  From  this,  we  con- 
clude that  the  convergence  constant  is  the  more  sensitive  input  parameter. 

Let  us  consider  another  second  order  dynamic  system  with  transfer  function, 


Hj{z)  = 


0.5-  0.2 1     +0.445  2 


-2 


1 


:"''  +  0.94z~2 


This  system  has  complex  poles  at  z  =  (0.5  +J0.8)  and  complex  zeros  located  at 
i  =  (0.2  ±j0.7).  Using  a  convergence  constant  p  =  0.005  and  power  level  weighting 
factor  p  =  0.97,  the  adaptive  ARMA  digital  lattice  filter  which  models  this  system  has 
steady-state  lattice  parameters  as  follows, 

terminal  condition  b0  =  0.500027 

lattice  coefficients  rl  =  0.104654 

r\  =  0.296658 

k\  =  -0.429232 

A?=  0.627718 
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Figure   16.      Second  order  ARMA  lattice  filter,  terminal  condition  unity. 
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Convergence  properties  of  this  adaptive  filter  are  shown  in  Figure  17.  In  this  example, 
convergence  was  obtained  after  approximately  2500  iterations.  Again,  by  changing  the 
values  of  u  and  p  slightly,  covergence  deteriorated  with  the  convergence  constant  /i  being 
the  more  sensitive  parameter. 

Next  we  consider  a  third  order  reference  system  with  known  transfer  function, 

0.5  -  0.95  z"1  +  1.33  z~2  -  0.979  z~3 


J'  1  -  1.69  z_1+ 0.962  z~2- 0.2  z~3 

The  adaptive  ARMA  digital  lattice  filter  which  describes  this  system  has  the  following 
steady-state  lattice  parameters, 

terminal  condition  b0  =  0.499970 

lattice  coefficients  r\  =  -0.328447 
r]  =  0.399472 
r\=  -0.652706 
k\  =  -O.S2173S 
k]=  -0.091111 
k\  =    —0.133333 

These  parameters  were  obtained  using  a  convergence  constant  n  =  0.015  and  power  level 
weighting  factor  p  =  0.9.  Convergence  properties  of  this  adaptive  filter  are  shown  in 
Figure  18.  Steady-state  values  for  the  lattice  coefficients  were  obtained  after  approxi- 
mately 7100  iterations.  It  is  reasonable  to  assume  that  a  third  order  system  will  converge 
more  slowly  than  a  second  order  system.  The  number  of  iterations  required  for  this  third 
order  system  to  converge  is  consistant  with  convergence  rates  of  other  adaptive  algo- 
rithms which  model  third  order  systems  [Ref.  16].  The  input  parameter  p  was  again 
found  to  be  the  more  sensitive  parameter. 

In  all  the  previous  examples,  the  values  of  p  and  p  may  or  may  not  be  optimum 
values.  That  is,  an  exhaustive  search  of  all  combinations  of  /i  and  p  was  not  performed 
to  demonstrate  convergence  of  the  algorithm.  Nevertheless,  a  number  of  different  ways 
of  realizing  the  value  of  the  convergence  constant  \x  have  been  reported  in  the  literature. 
In  one  method,  Mikhael  et.  al.  [Ref.  17]  have  obtained  a  variable  p  by  using  a  self  opti- 
mizing technique.  In  this  method,  p  is  calculated  from  the  input  data  as  an  iteration 
process  and  is  individually  determined  for  each  filter  parameter.  In  another  method  p  is 
chosen  by  using  a  variable  step  LMS  technique  [Ref.   18],  where  the  range  of  p  is 
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Figure   17.      Second  order  ARMA  lattice  filter,  terminal  condition  of  0.5. 
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Figure   18.      Third  order  ARMA  lattice  filter,  terminal  condition  of  0.5. 
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specified  by  p.m3X  and  ^min  which  are  within  the  bounds  described  by  equation  (3.4)  of 
Chapter  III.  These  techniques  of  choosing  /x  during  the  adaptation  process  have  been 
shown  to  improve  filter  convergence.  They,  however,  require  additional  computations 
to  achieve  this  faster  convergence.  When  a  combination  of  the  parameters  li  and  p  was 
obtained  which  yielded  convergence,  these  values  were  chosen  for  examples.  Besides  the 
examples  reported,  simulation  studies  have  been  carried  out  for  several  other  cases.  In 
all  cases,  however,  definite  convergence  of  the  algorithm  has  been  observed. 

In  summary,  we  have  demonstrated  through  computer  simulation  that  the  derived 
adaptive  lattice  algorithm  is  suited  for  system  identification  modeling.  Furthermore,  we 
have  shown  that  there  is  flexibility  in  choosing  the  values  of  the  convergence  constant 
ll  and  weighting  factor  p  .  Some  techniques  for  selecting  (computing)  the  value  of  the 
convergence  constant  have  been  introduced.  These  methods  improve  convergence  at  the 
cost  of  additional  computations. 

A.     CONCLUSIONS 

In  this  thesis  we  have  demonstrated  that  the  ARMA  parameter  estimation  algo- 
rithm proposed  in  [Ref.  4:  pp.  619-621]  is  a  valid  method  for  obtaining  approximations 
to  reference  models.  Furthermore,  the  criterion  used  to  derive  the  algorithm  is  a  gener- 
alized form  of  the  Mullis-Roberts  criterion  for  least  squares  modeling.  The  AR  and  MA 
parameters  of  the  ARMA  model  can  be  updated  independently  as  their  respective  orders 
increase  by  one.  From  the  recursive  prediction  error  formulas,  an  ARMA  digital  lattice 
filter  was  designed  with  arbitray  AR  and  MA  orders. 

For  the  ARMA  digital  lattice  filter,  we  derived  an  adaptive  lattice  algorithm.  This 
algorithm  was  based  on  the  least  mean  square  method  of  optimizing  coefficients.  The 
derived  adaptive  lattice  algorithm  can  easily  compute  the  values  of  the  lattice  coefficients 
from  available  data.  The  algorithm  simplified  the  number  of  coeflicents  required  to  be 
updated  from  eight  coefficients  per  elementary  lattice  section  to  only  two  such  that  the 
filter  can  be  completely  described.  This  savings  in  computational  effort  makes  the  al- 
gorithm attractive  for  identification  of  unknown  systems  since  many  systems  require  an 
ARMA  model  for  parsimonious  modeling. 

Convergence  of  the  adaptive  lattice  algorithm  was  demonstrated  with  several  exam- 
ples in  Chapter  III.  The  number  of  iterations  required  before  convergence  varied  greatly 
between  second  and  third  order  models  as  well  as  within  second  order  models.  Optimum 
convergence  rates  were  not  sought  after  as  much  as  proving  the  convergence  of  the  al- 
gorithm.   Rapid  convergence  rates  were  demonstrated  in  Chapter  II  for  a  second  order 
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system  upon  completion  of  an  extensive  search  for  the  optimun  values  of  the  conver- 
gence constant  and  power  level  weighting  factor. 

Although  a  method  of  converting  between  direct  form  and  lattice  realizations  for 
second  order  ARM  A  filters  was  developed,  a  general  algorithm  to  perform  this  trans- 
formation given  the  ARMA  lattice  filter  design  was  not  obtained.  When  solving  for  a 
transformation  between  filter  realizations  of  third  order  the  solution  is  hindered  by 
nonlinearities. 

The  objectives  of  the  thesis  were  sucsessfully  accomplished.  Some  suggestions  for 
future  work  include  the  following:  (i)  extensive  theoretical  analysis  for  determining  op- 
timum values  for  \i  ,  (ii)  derivation  of  a  generalized  algorithm  which  converts  any  given 
ARMA  transfer  function  into  a  set  of  lattice  parameters,  (Hi)  development  of  theoretical 
convergence  models  for  the  ARMA  adaptive  lattice  algorithm  and  analysis  of  these 
models  and  (iv)  application  of  the  adaptive  lattice  filter,  both  analysis  and  synthesis 
forms,  in  modeling  such  practical  signals  as  speech.  ARMA  lattice  filter  modeling  has 
considerable  application  potential  because  of  its  very  accurate  modeling  of  nearly  any 
sisnal  or  svstem  of  interest. 
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APPENDIX  A.     MAIN  PROGRAM  TO  ESTIMATE  ARMA  PARAMETERS 

C  THIS  PROGRAM  COMPUTES  THE  ARMA  PARAMETERS  AS  THE  AR  OR  MA  ORDER 

C  OF  AN  ARMA  MODEL  INCREASES  BY  ONE.  IT  USES  THE  ARMA  PARAMETER 

C  ESTIMATION  ALGORITHM  PROPOSED  BY  MIYANAGA,  NAGAI  AND  MIKI. 

C 

C  VARIABLE  DEFINITIONS 

C 

C  VN    -  INPUT  VECTOR  CONTAINING  DATA  GENERATED  AS  WHITE  GAUSSIAN 

C  NOISE 

C  Y       OUPUT  VECTOR  OF  DIFFERENCE  EQUATION  WITH  VN  AS  INPUT 

C  RX    -  AUTOCORRELATION  DATA  OF  INPUT  VN 

C  RY      AUTOCORRELATION  DATA  OF  OUTPUT  Y 

C  RXY   -  CROSS CORRELATION  DATA  OF  INPUT  AND  OUTPUT 

C  RYX    -  CROSSCORRELATION  DATA  OF  OUTPUT  AND  INPUT 

C  NDATA  -  NUMBER  OF  INPUT  DATA  POINTS 

C  KDATA  -  NUMBER  OF  INITIAL  DATA  POINTS  TO  DISREGARD 

C  XI     -  INPUT  DATA  VECTOR  AFTER  DISCARDING  KDATA  POINTS 

C  Yl    -  OUTPUT  DATA  VECTOR  AFTER  DISCARDING  KDATA  POINTS 

C  A     -  VECTOR  CONTAINING  CALCULATED  AR  PARAMETERS 

C  B     -  VECTOR  COINTAINING  CALCULATED  MA  PARAMETERS 

C  TXA    -  VECTOR  CONTAINING  COEFFICIENTS  FOR  FORWARD  PREDICTION  OF 

C  INPUT 

C  TYA   -  VECTOR  CONTAINING  COEFFICIENTS  FOR  FORWARD  PREDICTION  OF 

C  OUPUT 

C  TXB    -  VECTOR  CONTAINING  COEFFICIENTS  FOR  FORWARD  PREDICTION  OF 

C  INPUT 

C  TYB    -  VECTOR  CONTAINING  COEFFICIENTS  FOR  FORWARD  PREDICTION  OF 

C  OUTPUT 

C  GA      VECTOR  CONTAINING  COEFFICIENTS  FOR  BACKWARD  PREDICTION  OF 

C  INPUT 

C  GB    -  VECTOR  CONTAINING  COEFFICIENTS  FOR  BACKWARD  PREDICTION  OF 

C  INPUT 

C  ZTA    -  VECTOR  CONTAINING  AR  COEFFICIENTS  FOR  BACKWARD  PREDICTION 

C  OF  OUTPUT  VALUES 

C  ZTB    -  VECTOR  CONTAINING  MA  COEFFICIENTS  FOR  BACKWARD  PREDICTION 

C  OF  OUTPUT  VALUES. 

C  NI    -  LENGTH  OF  DATA  VECTOR  XI 

C  NK   -  LENGTH  OF  DATA  VECTOR  XI  MINUS  ONE  USED  TO  START 

C  CORRELATION  COMPUTATIONS. 

C  NS   -  DESIRED  AR  ORDER  OF  ARMA  MODEL. 

C  NT   -  DESIRED  MA  ORDER  OF  ARMA  MODEL 

C  KS    -  CURRENT  AR  ORDER  OF  UPDATE 

C  KT    -  CURRENT  MA  ORDER  OF  UPDATE 

C  VX   -  EXPECTED  VALUE  OF  PREDICTION  ERROR  FOR  INPUT  SQUARED. 

C  VY   -  EXPECTED  VALUE  OF  PREDICTION  ERROR  FOR  OUTPUT  SQUARED. 

C  VXY   -  EXPECTED  VALUE  OF  PRODUCT  OF  PREDICTION  ERROR  FOR  INPUT 
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C  AND  OUTPUT. 

C      VG     EXPECTED  VALUE  OF  BACKWARD  PREDICTION  ERROR  OF  INPUT 

C  SQUARED. 

C      VZ    -  EXPECTED  VALUE  OF  BACKWARD  PREDICTION  ERROR  OF  OUTPUT 

C  SQUARED. 

C      VGZ   -  EXPECTED  VALUE  OF  PRODUCT  BETWEEN  BACKWARD  PREDICTION 
ERRORS  OF  INPUT  AND  OUTPUT. 

C 

DIMENSION  VN(-5:  2006), Y( -5:  2006), RX(0:  2000), RY(0:  2000) ,RXY( 0:  2000) 
DIMENSION  RYX(0: 2000 ),X1( 2000), Yl( 2000), X( 12), A(0:  21),B(0:  21) 
DIMENSION  TXA(0:  21) ,TYA(0:  21) ,TXB(0:  21) ,TYB(0:  21) ,GA(0:  21) 
DIMENSION  GB(0:21),ZTA(0:  21),ZTB(0:  21) 

C 

C     INPUT  DATA  INFORMATION 

C 

WRITE  (6,1) 

1  FORMAT  (/'  ENTER  THE  NUMBER  OF  DATA  POINTS:  ') 
READ  (6,*)  NDATA 

WRITE  (6,2) 

2  FORMAT  (/'  ENTER  NUMBER  OF  INITIAL  DATA  POINTS  TO  DISREGARD: ') 
READ  (6,*)  KDATA 

C 

C 

C     INITIALIZE  ARRAYS 

C 

DO  10  L=-5, NDATA 

VN(L)=0 

Y(L)  =0 
10    CONTINUE 

DO  15  L=0,20 

A(L)=0 

B(L)=0 

TXA(L)=0 

TXB(L)=0 

TYA(L)=0 

TYB(L)=0 

GA(L)  =0 

GB(L)  =0 

ZTA(L)=0 

ZTB(L)=0 
15    CONTINUE 

DO  20  L=0,1999 

RX(L)=0 

RY(L)=0 

RXY(L)=0 

RYX(L)=0 
20    CONTINUE 

DO  25  L=l,12 
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X(L)=0 
25    CONTINUE 
C 

c 

C     GENERATE  WHITE  GAUSSIAN  NOISE  INPUT 
C 

ISIZE=NDATA 

IX=152255 

ISORT=0 

MUL=2 

DO  30  K=1,ISIZE 

CALL  SRND(IX,X,12,MUL,ISORT) 
XT=-6.  0 

DO  35  1=1,12 
35        XT=XT+X(I) 

VN(K)=XT 
30    CONTINUE 
C 

C 

C  COMPUTE  OUTPUT  OF  REFERENCE  MODEL  FILTER  AND  DISREGARD  SPECIFIED 

C  NUMBER  OF  DATA  POINTS 

C 

DO  40  L=1,NDATA 

C  Y(L)=VN(L)  +  1.  6*Y(L-l)-0.  95*Y(L-2) 

C  Y(L)=VN(L)+0.  2*Y(L-l)-0.  62*Y(L-2)+0.  152*Y( L-3) -0.  3016*Y(L-4) 

C  Y(L)=VN(L)-0.  2*VN(L-l)+0.  62*VN(L-2) -0.  152*VN(L-3)+0.  3016*VN(L-4) 

C  Y(L)=VN(L)+0.  2*VN(L-l)-0.  99*VN(L-2)+0.  2*Y(L-l)-0.  62*Y(L-2)+0.  152*Y 

C  &(L-3)-0.  3016*Y(L-4) 

C  Y(L)=VN(L)-1.6*VN(L-1)+1.  45*VN(L-2)+l.  2*Y(L-l)-0.  72*Y(L-2) 

C  Y(L)=VN(L)+0.  2*VN(L-l)-0.  35*VN(L-2)  +  l.  4*Y(  L-l)  -0.  85*Y(L-2) 

Y(L)=0.5*VN(L)-0.  2*VN(L-l)+0.445*VN(L-2)+Y(L-l)-0.  94*Y(L-2) 

C  Y(L)=VN(L)-2.  7*VN(L-l)+3.21*VN(L-2)-1.595*VN(L-3)  +  1.95*Y(L-l)-l.  62 

C  &*Y(L-2)+0. 54*Y(L-3) 

C  Y(L)=VN(L)-1.  0*VN(L-l)+0.  89*VN(L-2)+0.  40*Y( L- 1) -0.  2121*Y(L-2)-0.  20 

C  &894*Y(L-3)-l.  810373*Y(L-4) 

C  Y(L)=0.5*VN(L)-0.  95*VN(L-1)+1.  33*VN(L-2)-0.  979*VN(L-3)  +  l.  69*Y(L-1) 

C  &-0. 962*Y(L-2)+0.  20*Y(L-3) 

C  Y(L)=0.5*VN(L)-0.4*VN(L-l)+0.89*VN(L-2)  +  1.69*Y(L-l)-0.962*Y(L-2)+0 

C  &. 2*Y(L-3) 

C  Y(L)=0.5*VN(L)-0.4*VN(L-l)+0.89nW(L-2)+0.2*Y(L-l)+0.25*Y(L-2)-0.0 

C  &5*Y(L-3) 

C  Y(L)=0.  5*VN(L)-0.4*VN(L-l)+0.89*VN(L-2)+0.89*Y(L-l)-0.  25*Y(L-2) 

C  Y(L)=.  0154*VN(L)+.  0642*VN(L-l)+0.  0642*VN(L-2)+0.  0154*VN(L-3)+l.  99* 

C  &Y(L-1)-1.  57*Y(L-2)+0.  4583*Y(L-3) 

C  Y(L)=0.  5*VN(L)+0.  256*VN(L-l)+0.  1234*VN(L-2)+0.  0987*VN(L-3) 

40  CONTINUE 


52 


LJ=KDATA+1 

DO  45  L=LJ,NDATA 

LK=L-KDATA 

X1(LK)=VN(L) 

Y1(LK)=Y(L) 

45  CONTINUE 

C     WRITE  (*,77)  (Y1(K),  K=200,1800) 

77    FORMAT  (5( 1X,F10. 6) ) 

C 

C 

C     COMPUTE  AUTO -CORRELATION  AND  CROSS -CORRELATION  TERMS 

C 

46  NI=NDATA-KDATA 
NK=NI-1 

CALL  CORREL  (NI ,50 ,X1 , Yl ,RX,RY,RXY,RYX,NK) 

47  DO  50  L=0,10 

WRITE  (*,200)  RX(L),RY(L),RXY(L),RYX(L) 
WRITE  (9,200)  RX(L),RY(L),RXY(L),RYX(L) 
WRITE  (9,201) 
50    CONTINUE 

WRITE  (9,201) 

200  FORMAT  (2X,4(2X,F14. 9)) 

201  FORMAT  ('  ') 
C 

C 

C     INPUT  THE  DESIRED  AR  AND  MA  ORDERS  THEN  DEFINE  INITIAL  CONDITIONS 
C 

48  WRITE  (6,3) 

3  FORMAT  (/'  ENTER  THE  DESIRED  AR  ORDER:  ') 
READ  (6,*)  NS 

WRITE  (6,4) 

4  FORMAT  (/'  ENTER  THE  DESIRED  MA  ORDER: ') 
READ  (6,*)  NT 

C 

KS=0 

KT=0 

A(0)  =  1.0 

VX=RX(0) 

VY=RY(0) 

VXY=-RYX(0) 

VG=RX(0) 

VZ=RY(0) 

VGZ=-RYX(0) 

TXB(0)=1.0 

TYA(0)=1.0 

GB(0)  =1.  0 

ZTA(0)  =  1.0 
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c 

C     ESTIMATE  THE  ARMA  PARAMETERS 
C 

300  IF  (NT.  EQ.  0.  AND.KS.LT.  NS)  THEN 

KS=KS+1 
CALL  NEWAR(KS,KT,VX,VY,VXY,VG,VZ,VGZ,TXA,TXB,TYA,TYB,GA,GB,ZTA,ZTB 
&,RX,RY,RXY,RYX,A,B) 
GOTO  300 
ELSE 

301  IF  (NS.EQ.  0.  AND.KT.  LT.  NT)  THEN 

KT=KT+1 
CALL  NEWMA(KS,KT,VX,VY,VXY,VG,VZ,VGZ,TXA,TXB,TYA,TYB,GA,GB,ZTA,ZTB 
fie.RXjRY.RXY.RYX^B) 
GOTO  301 
ENDIF 
ENDIF 
C 

IF  (NS.NE.  O.OR.  NT.  NE.  0)  THEN 

302  IF  (NS.GE.NT.  AND.  NT.  NE.  0.  AND.  KT.  LT.  NT)  THEN 

KS=KS+1 
CALL  NEWAR( KS , KT , VX , VY , VXY , VG , VZ , VGZ , TXA , TXB , TYA , TYB , GA , GB , ZTA , ZTB 
&,RX,RY,RXY,RYX,A,B) 
KT=KT+1 
CALL  NEWMA( KS , KT , VX , VY , VXY , VG , VZ , VGZ , TXA , TXB , TYA , TYB , GA , GB , ZTA , ZTB 
&,RX,RY,RXY,RYX,A,B) 
GOTO  302 
ELSE 

303  IF  (KS.LT.  NS)  THEN 

KS=KS+1 
CALL  NEWAR(KS,KT,VX,VY,VXY,VG,VZ,VGZ,TXA,TXB,TYA,TYB,GA,GB,ZTA,ZTB 
&,RX,RY,RXY,RYX,A,B) 
GOTO  303 
ENDIF 
ENDIF 
ENDIF 
C 

C 

C     PRINT  ESTIMATED  ARMA  PARAMETERS 

C 

WRITE  (*,211) 

WRITE  (9,211) 

WRITE  (*,210)  (A(K),  K=1,KS) 

WRITE  (9,210)  (A(K),  K=1,KS) 

WRITE  (-,211) 

WRITE  (9,211) 
211   F0RMATC  ') 

WRITE  (*,210)  (B(K),  K=0,KT) 
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WRITE  (9,210)  (B(K),  K=0,KT) 
210   FORMAT  ( '  ' , IX, 4( 2X,F13. 10) ) 

STOP 

END 
C 

c 

C     SUBROUTINE  TO  COMPUTE  CORRELATION  TERMS 
C 

SUBROUTINE  C0RREL(N,LAG,X,Y,RX,RY,RXY,RYX,NK1) 
REAL  X(0:  NK1) ,Y(0:  NK1) ,RX(0:  2000) ,RY(0: 2000) ,RXY(0: 2000) 
REAL  RYX( 0: 2000 ) , SUM1 , SUM2 , SUM3 , SUM4 
DO  70  K=0,LAG 
NJ=N-1-K 
SUM1=0 
SUM2=0 
SUM3=0 
SUM4=0 
ANK=NJ 

DO  60  J=0,NJ 

SUM1=SUM1+X(J+K)*X(J) 
SUM2=SUM2+Y( J+K)*Y( J) 
SUM3=SUM3+X(  J+K)*Y( J) 
SUM4=SUM4+X(J)*Y(J+K) 
60        CONTINUE 

RX(K)=SUM1/ANK 
RY(K)=SUM2/ANK 
RYX(K)=SUM3/ANK 
RXY(K)=SUM4/ANK 
70    CONTINUE 
RETURN 
END 
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APPENDIX  B.     SUBROUTINE  FOR  MAIN  PROGRAM 


SUBROUTINE  NEWAR( KS , KT , VX , VY , VXY , VG , VZ , VGZ , TXA , TXB , TYA , TYB , GA , GB , Z 
&TA,ZTB,RX)RY,RXY,RYX,A,B) 

THIS  SUBROUTINE  COMPUTES  AR  PARAMETER  VALUES  FOR  AN  ARMA  MODEL  AS 
THE  AR  ORDER  INCREASES  BY  ONE. 

VARIABLE  DEFINITIONS 


CURRENT  AR  ORDER 

CURRENT  MA  ORDER 

AUTOCORRELATION  DATA  OF  INPUT  VN 

AUTOCORRELATION  DATA  OF  OUTPUT  Y 

CROSSCORRELATION  DATA  OF  INPUT  AND  OUTPUT 

CROSSCORRELATION  DATA  OF  OUTPUT  AND  INPUT 

VECTOR  CONTAINING  AR  COEFFICIENTS  FOR  FORWARD  PREDICTION 

OF  INPUT. 

VECTOR  CONTAINING  MA  COEFFICIENTS  FOR  FORWARD  PREDICTION 

OF  INPUT 

VECTOR  CONTAINING  AR  COEFFICIENTS  FOR  FORWARD  PREDICTION 

OF  OUTPUT 

VECTOR  CONTAINING  MA  COEFFICIENTS  FOR  FORWARD  PREDICTION 

OF  OUTPUT 

VECTOR  CONTAINING  AR  COEFFICIENTS  FOR  BACKWARD  PREDICTION 

OF  INPUT 

VECTOR  CONTAINING  MA  COEFFICIENTS  FOR  BACKWARD  PREDICTION 

OF  INPUT 

VECTOR  CONTAINING  AR  COEFFICIENTS  FOR  BACKWARD  PREDICTION 

OF  OUTPUT 

VECTOR  CONTAINING  MA  COEFFICIENTS  FOR  BACKYARD  PREDICTION 

OF  OUTPUT 

ARRAY  WHICH  STORES  CURRENT  VALUES  OF  TXA 

ARRAY  WHICH  STORES  CURRENT  VALUES  OF  TYA 

ARRAY  WHICH  STORES  CURRENT  VALUES  OF  GA 

ARRAY  WHICH  STORES  CURRENT  VALUES  OF  ZTA 

ARRAY  WHICH  STORES  CURRENT  VALUES  OF  TXB 

ARRAY  WHICH  STORES  CURRENT  VALUES  OF  TYB 

ARRAY  WHICH  STORES  CURRENT  VALUES  OF  GB 

ARRAY  WHICH  STORES  CURRENT  VALUES  OF  ZTB 

VECTOR  UPDATED  AR  COEFFICIENTS  OF  ARMA  MODEL 

VECTOR  CONTAINING  UPDATED  MA  COEFFICIENTS  OF  ARMA  MODEL. 

CONSTANT  COMPUTED  FROM  CORRELATION  DATA  AND  PREDICTION 

ERROR  COEFFICIENTS. 

CONSTANT  COMPUTED  FROM  CORRELATION  DATA  AND  PREDICTION 

ERROR  COEFFICIENTS. 

CONSTANT  COMPUTED  FROM  CORRELATION  DATA  AND  PREDICTION 


c 

c 

THIS 

c 

c 

c 

c 

c 

KS 

c 

KT 

c 

RX 

c 

RY 

c 

RXY 

c 

RXY 

c 

TXA 

c 

c 

TXB 

c 

c 

TYA 

c 

c 

TYB 

c 

c 

GA 

c 

c 

GB 

c 

c 

ZTA 

c 

c 

ZTB 

c 

c 

C 

c 

D 

c 

E 

c 

F 

c 

P 

c 

Q 

c 

R 

c 

S 

c 

A 

c 

B 

c 

TAU1 

c 

c 

TAU2 

c 

c 

TAU3 
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c 

c 

TAU4 

c 

c 

XMU1 

c 

YMU1 

c 

MU2 

c 

XMU3 

c 

YMU3 

c 

MU4 

c 

MU5 

c 

DET 

c 

c 

ERR 

c 

10 
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ERROR  COEFFICIENTS. 
CONSTANT  COMPUTED  FROM 
ERROR  COEFFICIENTS. 
COEFFICIENT  OF  AR-TYPE 
COEFFICIENT  OF  AR-TYPE 
COEFFICIENT  OF  AR-TYPE 
COEFFICIENT  OF  AR-TYPE 
COEFFICIENT  OF  AR-TYPE 
COEFFICIENT  OF  AR-TYPE 
COEFFICIENT  OF  AR-TYPE 
DETERMINANT  OF  PREDICTI 
VX,VY,VXY. 

ERROR  BETWEEN  REFERENCE 
REALIZATION  OUTPUT. 


CORRELATION  DATA  AND  PREDICTION 

RECURSIVE  FORMULA 
RECURSIVE  FORMULA 
RECURSIVE  FORMULA 
RECURSIVE  FORMULA 
RECURSIVE  FORMULA 
RECURSIVE  FORMULA 
RECURSIVE  FORMULA 
ON  ERROR  MATRIX  COMPOSED  OF 

MODEL  OUTPUT  AND  LATTICE 


DIMENSION  RX(0: 2000), RY(0:  2000),RXY(0:  2000),RYX(0:  2000), A(0:  21) 
DIMENSION  B(0: 21) ,TXA( 0: 21) ,TYA( 0:  21) ,TXB(0:  21) ,TYB( 0:  21),GA(0:21) 
DIMENSION  GB(0:21),ZTA(0:  21),ZTB(0:  21)  ,C(0:  21)  ,D(0:  21)  ,E(  0:  21) 
DIMENSION  P(0:  21),Q(0:  21),R(0:  21),S(0:  21),F(0:21) 
REAL  MU5,MU2,MU4 

COMPUTE  VALUES  FOR  TAU1  THROUGH  TAU4 

T1S=0 

T2S=0 

T3S=0 

T4S=0 

KI=KS-1 

DO  10  1=0, KI 

T1S=T1S-RYX(I+1)*ZTA(KI-I) 

T2S=T2S+RY(I+1)*ZTA(KI-I) 

T3S=T3S+RY(I+1)*GA(I) 

T4S=T4S+RY(I+1)*ZTA(I) 
CONTINUE 
T1T=0 
T2T=0 
T3T=0 
T4T=0 
DO  20  J=0,KT 

T1T=T1T+RX(J+1)*ZTB(KT-J) 

T2T=T2T-RXY(  J+l )*ZTB( KT- J) 

T3T=T3T-RYX(KI-KT+1+J)*GB(J) 

T4T=T4T-RYX( KI -KT+1+J)*ZTB( J) 
CONTINUE 
TAU1=T1S+T1T 
TAU2=T2S+T2T 
TAU3=T3S+T3T 
TAU4=T4S+T4T 
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C     COMPUTE  VALUES  FOR  THE  REFLECTION  COEFFICIENTS 
C 

XMU1=-TAU1/VZ 

YMU1=-TAU2/VZ 

MU2  =-VGZ/VZ 

DET  =VX*VY-VXY*VXY 

XMU3= - ( VY*TAU 1 - VXY*TAU2 ) /DET 

YMU3=-(-VXY*TAUl+VX*TAU2)/DET 

MU4  =( VGZ*TAU4 - VZ*TAU3 ) / ( VGZ*VGZ-VG*VZ ) 

MU5  =MU4*MU2 
C 
C     COMPUTE  ARMA  COEFFICIENTS 

c 

DO  16  K=0,KS 

C(K)=TXA(K) 

D(K)=TYA(K) 

E(K)=GA(K) 

F(K)=ZTA(K) 
16    CONTINUE 

DO  45  J=1,KS 
TXA(J)=C(J)+XMU1*F(KS-J) 
TYA(J)=D(J)+YMU1*F(KS-J) 
GA(J)=E(J-1)+MU2*F(J-1) 

ZTA(J)=F(J)+XMU3*C(KS-J)+YMU3*D(KS-J)+MU4*E(J-1)+MU5*F(J-1) 
45    CONTINUE 

DO  31  K=0,KT 

P(K)=TXB(K) 

Q(K)=TYB(K) 

S(K)=ZTB(K) 

R(K)=GB(K) 
31    CONTINUE 

DO  55  J=1,KT 
TXB(J)=P(J)+XMU1*S(KT+1-J) 
TYB(J)=Q(J)+YMU1*S(KT+1-J) 
GB(J)  =R(J)+MU2*S(J) 

ZTB(J)=S(J+1)+XMU3*P(KT-J)+YMU3*Q(KT-J)+MU4*R(J)+MU5*S(J) 
55    CONTINUE 
C     WRITE  (",176)  KS 
C     WRITE  (9,176)  KS 
176   FORMAT(I2) 

C     WRITE  (*,175)  (ZTB(K),  K=0,KT) 
C     WRITE  (9,175)  (ZTB(K),  K=0,KT) 
175   FORMAT  (4( 1X,F10. 5) ) 
C 

C     UPDATE  ERRORS 
C 
650   FORMAT  (/'  S  UPDATE  ERROR  IS:  ' ,F15. 10) 

VX  =VX+XMU1*TAU1 

VY  =VY+YMU1*TAU2 

VXY=VXY+XMU1*TAU2 
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VG  =VG+MU2*VGZ 

VZ=VZ+(XMU3*TAU1+YMU3*TAU2)+MU4*TAU3+MU5*TAU4 

VGZ=TAU3+MU2*TAU4 

ERR=VY-(VXY**2)/VX 

WRITE  (*,650)  ERR 

WRITE  (9,650)  ERR 

WRITE  (*,888)  VX,VY,VG,VZ 
888   FORMAT  (4( 1X,F10.  6) ) 
C 

C     COMPUTE  MODEL  COEFFICIENTS 
C 

DO  65  J=1,KS 
A(J)=TYA(J)-TXA(J)*VXY/VX 
65    CONTINUE 

DO  70  J=1,KT 
B(J)=TYB(J)-TXB(J)*VXY/VX 
70    CONTINUE 

B(0)=-VXY/VX 

RETURN 

END 
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APPENDIX  C.     SUBROUTINE  FOR  MAIN  PROGRAM 


SUBROUTINE  NEWMA(KS ,KT, VX, VY , VXY, VG, VZ , VGZ , TXA, TXB , TYA, TYB ,GA,GB , Z 
&TA,ZTB,RX,RY,RXY,RYX,A,B) 

THIS  SUBROUTINE  COMPUTES  MA  PARAMETER  VALUES  FOR  AN  ARMA  MODEL  AS 
THE  MA  ORDER  INCREASES  BY  ONE. 

VARIABLE  DEFINITIONS 


CURRENT  AR  ORDER 

CURRENT  MA  ORDER 

AUTOCORRELATION  DATA  OF  INPUT  VN 

AUTOCORRELATION  DATA  OF  OUTPUT  Y 

CROSSCORRELATION  DATA  OF  INPUT  AND  OUTPUT 

CROSSCORRELATION  DATA  OF  OUTPUT  AND  INPUT 

VECTOR  CONTAINING  AR  COEFFICIENTS  FOR  FORWARD  PREDICTION 

OF  INPUT. 

VECTOR  CONTAINING  MA  COEFFICIENTS  FOR  FORWARD  PREDICTION 

OF  INPUT 

VECTOR  CONTAINING  AR  COEFFICIENTS  FOR  FORWARD  PREDICTION 

OF  OUTPUT 

VECTOR  CONTAINING  MA  COEFFICIENTS  FOR  FORWARD  PREDICTION 

OF  OUTPUT 

VECTOR  CONTAINING  AR  COEFFICIENTS  FOR  BACKWARD  PREDICTION 

OF  INPUT 

VECTOR  CONTAINING  MA  COEFFICIENTS  FOR  BACKWARD  PREDICTION 

OF  INPUT 

VECTOR  CONTAINING  AR  COEFFICIENTS  FOR  BACKWARD  PREDICTION 

OF  OUTPUT 

VECTOR  CONTAINING  MA  COEFFICIENTS  FOR  BACKWARD  PREDICTION 

OF  OUTPUT 

ARRAY  WHICH  STORES  CURRENT  VALUES  OF  TXA 

ARRAY  WHICH  STORES  CURRENT  VALUES  OF  TYA 

ARRAY  WHICH  STORES  CURRENT  VALUES  OF  GA 

ARRAY  WHICH  STORES  CURRENT  VALUES  OF  ZTA 

ARRAY  WHICH  STORES  CURRENT  VALUES  OF  TXB 

ARRAY  WHICH  STORES  CURRENT  VALUES  OF  TYB 

ARRAY  WHICH  STORES  CURRENT  VALUES  OF  GB 

ARRAY  WHICH  STORES  CURRENT  VALUES  OF  ZTB 

VECTOR  UPDATED  AR  COEFFICIENTS  OF  ARMA  MODEL 

VECTOR  CONTAINING  UPDATED  MA  COEFFICIENTS  OF  ARMA  MODEL. 

CONSTANT  COMPUTED  FROM  CORRELATION  DATA  AND  PREDICTION 

ERROR  COEFFICIENTS. 

CONSTANT  COMPUTED  FROM  CORRELATION  DATA  AND  PREDICTION 

ERROR  COEFFICIENTS. 

CONSTANT  COMPUTED  FROM  CORRELATION  DATA  AND  PREDICTION 


c 

c 

THIS 

c 

c 

c 

c 

c 

KS 

c 

KT 

c 

RX 

c 

RY 

c 

RXY 

c 

RXY 

c 

TXA 

c 

c 

TXB 

c 

c 

TYA 

c 

c 

TYB 

c 

c 

GA 

c 

c 

GB 

c 

c 

ZTA 

c 

c 

ZTB 

c 

c 

C 

c 

D 

c 

E 

c 

F 

c 

P 

c 

Q 

c 

R 

c 

S 

c 

A 

c 

B 

c 

TAU1P 

c 

c 

TAU2P 

c 

c 

TAU3P 
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c 

c 

TAU4P 

c 

c 

XETA1 

c 

YETA1 

c 

ETA2 

c 

XETA3 

c 

YETA3 

c 

ETA4 

c 

ETA5 

c 

DET 

c 

c 

ERR 

c 

ERROR  COEFFICIENTS. 

CONSTANT  COMPUTED  FROM  CORRELATION  DATA  AND  PREDICTION 

ERROR  COEFFICIENTS. 

COEFFICIENT  OF  MA-TYPE  RECURSIVE  FORMULA 

COEFFICIENT  OF  MA-TYPE  RECURSIVE  FORMULA 

COEFFICIENT  OF  MA-TYPE  RECURSIVE  FORMULA 

COEFFICIENT  OF  MA-TYPE  RECURSIVE  FORMULA 

COEFFICIENT  OF  MA-TYPE  RECURSIVE  FORMULA 

COEFFICIENT  OF  MA-TYPE  RECURSIVE  FORMULA 

COEFFICIENT  OF  MA-TYPE  RECURSIVE  FORMULA 

DETERMINANT  OF  PREDICTION  ERROR  MATRIX  COMPOSED  OF 

VX,VY,VXY. 

ERROR  BETWEEN  REFERENCE  MODEL  OUTPUT  AND  LATTICE 

DIMENSION  RX(0:2000),RY(0:  2000) ,RXY(0:  2000) ,RYX(0:  2000) ,A(0:  21) 
DIMENSION  B(0: 21),TXA(0: 21),TXB(0: 21),TYA(0: 21),TYB(0: 21),GA(0: 21) 
DIMENSION  GB(0: 21),ZTA(0: 21),ZTB(0:  21) ,C(0:  21) ,D( 0:  21) 
DIMENSION  E(0:  21),F(0:  21),P(0:  21),Q(0:  21)  ,R(0:  21)  ,S(0:  21) 
C 

C     COMPUTE  VALUES  FOR  TAU1  PRIME  THROUGH  TAU4  PRIME 
C 

T1TP=0 
T2TP=0 
T3TP=0 
T4TP=0 
KJ=KT-1 
DO  10  1=0, KJ 
T1TP=T1TP+RX(I+1)*GB(KJ-I) 
T2TP=T2TP-RXY(I+1)*GB(KJ-I) 
T3TP=T3TP+RX(I+1)*ZTB(I) 
T4TP=T4TP+RX ( I + 1 ) *GB ( I ) 
10    CONTINUE 
T1SP=0 
T2SP=0 
T3SP=0 
T4SP=0 

DO  20  J=0,KS 
T1SP=T1SP-RYX(J+1)*GA(KS-J) 
T2SP=T2SP+RY(J+1)*GA(KS-J) 
T3SP=T3SP-RXY(KJ-KS+1+J)*ZTA(J) 
T4SP=T4SP-RXY(KJ-KS+1+J)*GA(J) 
20    CONTINUE 

TAU1P=T1TP+T1SP 
TAU2P=T2TP+T2SP 
TAU3P=T3TP+T3SP 
TAU4P=T4TP+T4SP 
C 

C     COMPUTE  VALUES  FOR  REFLECTION  COEFFICIENTS 
C 

XETA1=-TAU1P/VG 
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YETA1=-TAU2P/VG 

ETA2  =-VGZ/VG 

DET  =VX*VY-VXY*VXY 

XETA3=-(VY*TAU1P-VXY*TAU2P)/DET 

YETA3=- ( -VXY*TAU1P+VX*TAU2P) /DET 

ETA4  =(VGZ*TAU4P-VG*TAU3P)/(VGZ*VGZ-VG*VZ) 

ETA5  =ETA4*ETA2 
C 

C     COMPUTE  ARMA  PARAMETERS 
C 

DO  16  K=0,KT 

P(K)=TXB(K) 

Q(K)=TYB(K) 

R(K)=GB(K) 

S(K)=ZTB(K) 
16    CONTINUE 
165   FORMAT  (4( 1X,F10. 5) ) 

DO  45  J=1,KT 
TXB(J)=P(J)+XETA1*R(KT-J) 
TYB(J)=Q(J)+YETA1*R(KT-J) 

GB(J)=R(J)+XETA3*P(KT-J)+YETA3*Q(KT-J)+ETA4*S(J-1)+ETA5*R(J-1) 
ZTB(J)=S(J-1)+ETA2*R(J-1) 
45    CONTINUE 

DO  41  K=0,KS 

C(K)=TXA(K) 

D(K)=TYA(K) 

E(K)=GA(K) 

F(K)=ZTA(K) 
41    CONTINUE 
175    FORMAT  (5(1X,F10. 5)) 

DO  55  J=1,KS 
TXA(J)=C(J)+XETA1*E(KS+1-J) 
TYA(J)=D(J)+YETA1*E(KS+1-J) 

GA(J)=E(J+1)+XETA3*C(KS-J)+YETA3*D(KS-J)+ETA4*F(J)+ETA5*E(J) 
ZTA(J)=F(J)+ETA2*E(J) 
55    CONTINUE 
C 

C     UPDATE  ERRORS 
C 

VX=VX+XETA1*TAU1P 

VY=VY+YETA1*TAU2P 

VXY=VXY+XETA 1*TAU2P 

VG=VG+(TAU1P*XETA3+TAU2P*YETA3)+ETA4*TAU3P+ETA5*TAU4P 

VZ=VZ+ETA2*VGZ 

VGZ=TAU3P+ETA2*TAU4P 

ERR=VY-(VXY**2)/VX 

WRITE  (*,66)  ERR 

WRITE  (9,66)  ERR 
66    FORMAT  (/*  T  UPDATE  ERROR  IS:  ' ,F15. 10) 
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WRITE  (*,889)  VX,VY,VG,VZ 
889   FORMAT  (4(1X,F10.6)) 
C 

C     COMPUTE  MODEL  COEFFICIENTS 
C 

DO  65  J=1,KT 
B(J)=TYB(J)-TXB(J)*VXY/VX 
65    CONTINUE 

DO  70  J=1,KS 
A(J)=TYA(J)-TXA(J)*VXY/VX 
70    CONTINUE 

B(0)=-VXY/VX 

RETURN 

END 
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APPENDIX  D.     ADAPTIVE  LATTICE  ALGORITHM  PROGRAM 

THIS  PROGRAM  CALCULATES  VALUES  OF  THE  LATTICE  COEFFI CENTS AND 
OUTPUT  OF  AN  ARMA  DIGITAL  LATTICE  FILTER  USING  AN  ADAPTIVE 

LATTICE  ALGORITHM. 

VARIABLE  DEFINITIONS 

X     -  ARRAY  OF  INPUT  DATA,  COMPUTER  GENERATED  WHITE  GAUSSIAN 

NOISE  WITH  UNIT  VARIANCE. 
AK    -  ARRAY  OF  LATTICE  COEFFICIENTS. 
R     -  ARRAY  OF  LATTICE  COEFFICIENTS. 
BO     -  TERMINAL  CONDITION  OF  LATTICE  REALIZATION. 
M     -  NUMBER  OF  LATTICE  STAGES  (EQUIVALENT  TO  ORDER  OF  ARMA 

MODEL). 
EXF    -  ARRAY  OF  FORWARD  PREDICTION  ERRORS  FOR  INPUT  X. 
EXB    -  ARRAY  OF  BACKWARD  PREDICTION  ERRORS  FOR  INPUT  X. 
EXBD   -  ARRAY  OF  DELAYED  BACKWARD  PREDICTION  ERRORS  FOR  INPUT  X. 
EYF    -  ARRAY  OF  FORWARD  PREDICTION  ERRORS  FOR  OUTPUT  Y. 
EYB   -  ARRAY  OF  BACKWARD  PREDICTION  ERRORS  FOR  OUTPUT  Y. 
EYBD   -  ARRAY  OF  DELAYED  BACKWARD  PREDICTION  ERRORS  FOR  OUTPUT  Y. 
ERROR  -  DIFFERENCE  BETWEEN  REFERENCE  MODEL  OUTPUT  AND  LATTICE 

REALIZATION  OUTPUT. 
YE    -  ARRAYS  CONTAINING  LATTICE  COEFFICIENT  VALUES  AT  EACH 

ITERATION. 
MU    -  CONVERGENCE  CONSTANT. 
RHO    -  WEIGHT  GIVEN  TO  CUURENT  POWER  LEVEL  AT  EACH  STAGE  OF  THE 

LATTICE  STRUCTURE. 
SIGK   -  POWER  LEVEL  USED  TO  NORMALIZE  CONVERGENCE  CONSTANT  WHEN 

UPDATING  AK  LATTICE  COEFFICIENTS. 
SIGR   -  POWER  LEVEL  USED  TO  NORMALIZE  CONVERGENCE  CONSTANT  WHEN 

UPDATING  R  LATTICE  COEFFICIENTS. 
SIGB   -  POWER  LEVEL  USED  TO  NORMALIZE  CONVERGENCE  CONSTANT  WHEN 

UPDATING  TERMINAL  CONDITION  BO. 

DIMENSION  EXF( 10) ,EXB( 10) ,EXBD( 10) ,EYF( 10) ,EYB( 10) ,EYBD( 10) ,R( 10) 

DIMENSION  AK(10),  X(9900),  ERROR(9900),  V(12),  YE(6,9900) 

REAL  MU 

SIGK=1. 

SIGR=1. 

SIGB=1. 

INITIALIZE  ARRAYS 

DO  5  1=1,10 

EXF(I)=0 

EXB(I)=0 


64 


EXBD(I)=0 

EYF(I)=0 

EYB(I)=0 

EYBD(I)=0 

R(I)=0 

AK(I)=0 

5  CONTINUE 

DO  6  1=1,9000 

X(I)=0 

ERROR(I)=0 

6  CONTINUE 
C 

C     ENTER  VALUE  OF  THE  CONVERGENCE  CONSTANT  MU  AND  VALUE  OF  RHO. 
C 

M=2 

N=300 

WRITE  (6,*)  'ENTER  MU1 

RE AD (6,*)  MU 

WRITE  (6,*)  'ENTER  RHO' 

READ  (6,*)  RHO 
C 

C     GENERATE  WHITE  GAUSSIAN  NOISE 
C 

ISIZE  =  N 

IX  =  152255 

I SORT  =  0 

MUL  =  2 

DO  7  K=  1,  ISIZE 
CALL  SRND(IX,V,12,MUL,IS0RT) 
XT=-6. 0 

DO  8  1=1,12 
8        XT=XT+V(I) 
X(K)=XT 

7  CONTINUE 
C 

C     COMPUTE  OUTPUT  OF  REFERENCE  MODEL  FILTER  AND  LATTICE  STRUCTURE 

C     THEN  COMPUTE  THE  ERROR. 

C 

C     REFERENCE  MODEL 

Y3=0 
Y2=0 
Y1=0 
X3=0 
X2=0 
X1=0 
B0=1. 

DO  100  1=1, N 
C     YF=X(I)-0.  8*X1+1.  78*X2+0.  89*Yl-0.  25*Y2 
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YF=0. 5*X(I)-0.  4*X1+.  89*X2+0.  89*Yl-0.  25*Y2 
C     YF=X(I)-2.  7*Xl+3.  21*X2-1.  595*X3+1.  95*Y1-1.  62*Y2+0.  54*Y3 
C     YF=0.  5*X(I)-0.  95*X1+1.  33*X2-0.  979*X3+1.  69*Yl-0.  962*Y2+0.  2*Y3 
C     YF=X(I)+0.  2*Xl-0.  35*X2+1.  4*Yl-0.  85*Y2 
C     YF=0.  5*X(I)-0.  2*Xl+0.  445*X2+1.  0*Yl-0.  94*Y2 

C     Y3=Y2 

Y2=Y1 

Y1=YF 
C     X3=X2 

X2=X1 

X1=X(I) 

c 

C     LATTICE  FILTER 
C 

EXF(1)=X(I) 

EXB(1)=X(I) 

DO  10  K=1,M 
10    EXF(K+1)=EXF(K)+R(K)*EXBD(K)-AK(K)*EYBD(K) 

EYF(M+1)=B0*EXF(M+1) 

DO  20  K=1,M 
20    EYF(M+l-K)=EYF(M+2-K)+R(M+l-K)*EXBD(M+l-K)-AK(M+l-K)*EYBD(M+l-K) 

EYB(1)=EYF(1) 

DO  30  K=1,M-1 

EXB(K+1)=EXBD(K)+R(K)*EXF(K)-R(K)*EYF(K) 
30    EYB(K+1)=EYBD(K)+AK(K)*EYF(K)-AK(K)*EXF(K) 

YL=EYF(1) 

ERROR(I)=YF-YL 

ERR=YF-YL 

CALL  UPDATE  (R,AKsEYBD,EXBD,ERR,MU,M,SIGK,SIGR,BO,RHO) 

CSB=EXF(M+1)*EXF(M+1) 

SIGB=RHO--'--SIGB+(  l-RHO)*CSB 

B0=B0+CMU/SIGB)*ERR*EXF(M+1) 

DO  40  K=1,M 

EXBD(K)=EXB(K) 
40    EYBD(K)=EYB(K) 

DO  50  J=1,M 

YE(J,I)=AK(J) 
50    YE(J+M,I)=R(J) 
202   FORMAT  (2(1X,F10. 6)) 
100   CONTINUE 
C 

C     PRINT  THE  ERROR  AND  VALUES  OF  THE  LATTICE  COEFFICIENTS. 
C 

WRITE  (''-,200)  (ERROR(K),  K=1,N,10) 

WRITE  (9,200)  (ERROR(K),  K=1,N,10) 
200   FORMAT  (5( 1X,F10. 6) ) 

WRITE  (9,209) 
209   FORMAT( '  ' ) 

WRITE  (*,201)  (R(K),  K=1,M) 
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WRITE  (*,201)  (AK(K),  K=1,M) 

WRITE  (9,201)  (R(K),  K=1,M) 

WRITE  (9,201)  (AK(K),  K=1,M) 

WRITE  (*,205)  BO 

WRITE  (9,205)  BO 
205    FORMAT  (F10. 6) 
201   FORMAT  (5(1X,F10. 6)) 
C 

C     CALL  PLOTTING  ROUTINES  TO  PLOT  ERROR  AND  LATTICE  COEFFICIENTS. 
C 

CALL  PLOT  ( ERROR, N) 
C     CALL  PLOT1  (YE,N) 

STOP 

END 
C 

C     SUBROUTINE  WHICH  UPDATES  LATTICE  COEFFICIENTS. 
C 

SUBROUTINE  UPDATE(R,AK,EYBD,EXBD,ERR,MU,M,SIGK,SIGR,B0,RHO) 

DIMENSION  R( 10) ,AK( 10) ,EYBD( 10) ,EXBD( 10) 

REAL  MU 

CSK=0. 

CSR=0. 

DO  20  J=1,M 

CSK=CSK+EYBD(J)*EYBD(J)*B0**2 
20    CSR=CSR+EXBD(J)*EXBD(J)*B0**2 

SIGK=RHO*SIGK+(l-RHO)*CSK 

SIGR=RHO*SIGR+( l-RHO)*CSR 

DO  10  J=1,M 

R(J)=R(J)+(MU/SIGR)*ERR*EXBD(J)*BO 

AK(J)=AK(J)-(MU/SIGK)*ERR*EYBD(J)*BO 
10    CONTINUE 

RETURN 

END 
C 

C     PLOTTING  ROUTINE  TO  PLOT  ERROR 
C 

SUBROUTINE  PLOT(Y,N) 

DIMENSION  Y(N),X(9900) 

DO  10  J=1,N 
10    X(J)=J 

CALL  TEK618 
C     CALL  PRTPLT(72,6) 
C     CALL  SHERPAC ADAPTIVE' ,' A' ,3) 

CALL  RESET('ALL') 

CALL  PAGE(  8.  50,6.0) 

CALL  HWROT('AUTO') 

CALL  XI NT AX 

CALL  AREA2D(5.0,3.  0) 

CALL  HEIGHT(0. 14) 

CALL  COMPLX 
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CALL  SHDCHR( 90.  0,1,0.  002,1) 

CALL  HEADIN(' LEARNING  CURVES$ ' , 100 , 2.  0 , 1) 

CALL  XNAMEC ' ITERATIONS$ ' , 100) 

CALL  YNAMEC 'ERROR $' ,100) 
C     CALL  MESSAGC   ADAPTIVE  FILTER   $',100,3.0,-0.8) 

CALL  THKFRM(0.  03) 

CALL  FRAME 

CALL  GRAF(0,' SCALE' ,N,-3.  00, 'SCALE'  ,3.  00) 
C     CALL  THKCRV(0.  02) 

CALL  CURVE(X,Y,N,0) 

CALL  ENDPL(O) 

CALL  DONEPL 

RETURN 

END 
C 

C     SUBROUTINE  TO  PLOT  LATTICE  COEFFICIENTS 
C 

SUBROUTINE  PL0T1  (YE,N) 

DIMENSION  YE(6,9900),X(9900),Y(9900),YD(9900),A(10) 

C TRUE  VALUES  OF  THE  PARAMETERS 

C     A(l)  =  -0.  240719 
C     A(2)=0.  125 
C     A(3)  =  -0.  195719 
C     A(4)=0.8900 
C     A(5)=0.S9 

DO  10  J=1,N 
10    X(J)=J 

CALL  TEK618 
C     CALL  PRTPLT(72,6) 
C     CALL  SHERPA('MENNECKE' ,'A' ,3) 
C PRINT  SHERPA  FILE:  SHERPA  XXYYZZXX  SHGRAPH  A 

CALL  RESET('ALL') 

CALL  PAGE(8.50,6.  0) 

CALL  HWROT( ' AUTO ' ) 

CALL  XINTAX 

CALL  AREA2D(5.  0,3.  0) 

CALL  HEIGHT(0.  14) 

CALL  COMPLX 

CALL  SHDCHR(90.  0,1,0.  002,1) 

CALL  HEADINC  PARAMETERS $'  ,100,2.0,1) 

CALL  XNAME(' ITERATIONS?' ,100) 

CALL  YNAME( 'MAGNITUDE?' ,100) 
C     CALL  MESSAG( 'ADAPTIVE  ARMA  LATTICE$ ' , 100,3.  0, -0.  8) 

CALL  THKFRM(0.  03) 

CALL  FRAME 

CALL  GRAF(0, 'SCALE' ,N,-1. 00, 'SCALE' ,1.  0) 
C     CALL  THKCRVC0.02) 
C TO  PLOT  ESTIMATES 

DO  20  K=l,4 
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DO  30  J=1,N 
30    Y(J)=YE(K,J) 

CALL  CURVE(X,Y,N,0) 
20    CONTINUE 
C TO  PLOT  TRUE  PARAMETERS 

DO  40  K=l,4 

DO  50  J=1,N 
50    Y(J)=A(K) 

CALL  DASH 

CALL  CURVE(X,Y,N,0) 
40    CONTINUE 

CALL  ENDPL(O) 

CALL  DONEPL 

RETURN 

END 
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APPENDIX  E.     DERIVATION  OF  DIFFERENCE  EQUATION  FOR  TWO 

STAGE  ARMA  LATTICE  FILTER 

The  lattice  filter  is  described  by  the  expressions  for  the  forward  and  backward  pre- 
diction errors,  equations  (3.6)  and  (3.7)  respectively, 


exh  (A)  =  x(k)  +  w\  x(k  -  1)  -  w\  y(k  -  1) 

exh  (k)  =  exh  (k)  +  w\  exbi  (k  -  1)  -  w\  ^  (k  -  I) 

e£  (A)  =  x(k  -  1)  +  w\  x(k)  -  w\  y{k)  (E-  1) 

eii(k)=y(k-\)-»-lx(k)  +  vl}ik) 

$  (k)  =  exh  (A)  +  ;v42  e*  (k  -  1)  -  w\  ^  (k  -  1) 

and  the  expression  for  the  filter  output, 

y(k)  =  e]^  (k)  +  w\  x(k  -  1 )  -  wl  y(k  -  1 )  (£  -  2) 

Substituting  for  ej  (A)  in  equation  (E-2)  yields, 

y(k)  =  ef2  (k)  +  vv42  «*  (A  -  1)  -  w]  «J  (A-  -  1)  +  w\  x(k  -  1)  -  iirj  >-(A  -  1)     (£-3) 

Substituting  for  ef  (k  —  1)  and  eL  {k  —  1)  in  equation  (E-3) 

y(k)  =  exh  (k)  +  u-2  [  x(k  -  2)  +  wj  *(A  -  1)  -  w\  y(k  -  1)  ] 

-  w2  [  j-(A  -  2)  -  w\  x(k-l)  +  wl  y(k-  1)  ]  (£  -  4) 

+  11.4  .x(A  —  1)  —  vv3  y(A  —  1) 

Substituting  for  e)  [k)  in  (E-4)  we  obtain 

y(k)  =  exh  (k)  +  wl  ex:  (k  -  1)  -  w2  e£  (A  -  1)  +  w42  *(A  -  2)  +  w\  vt-s1  x(k  -  1) 

— wA  wl  y(k  —  1)  —  w3  y(k  —  2)  +  w3  w^  x(k  —  1)  —  Uj  w2  y(k  —  1)  (£  —  5) 

+wlx(k-\)-W]3y(k-\) 

From  (E-l),  substitute  for  e$  (A  -  1)  and  e%  (A  -  1)  in  (E-5)  to  obtain 
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y{k)  =  exh  (A)  4  wj  [  x(k  -  2)  4-  w]  x(k  -  1)  -  wj  y(k  -  1)  ] 
-  w\  [  y(k  -  2)  -  w\  x{k  -1)4  w]  y{k  -  1)  ] 

4-  w4  x(A  -  2)  4-  w4  w2  x(k  -  1)  -  vv4  wj  v(A  -  1)  (£  -  6) 

-w3  >•(*  -  2)  +  w3  wj1  .v(A  -  1)  -  w3  wj  j-(A  -  1) 
4-  w4  .x(A'  —  1)  —  vv3  j(A  —  1) 

Now,  from  (E-l),  substituting  for  e)  (k)  in  (E-6)  yields 

y{k)  =  x{k)  4-  w\  x{k  -  1)  -  w\  y{k  -  1)  +  w2  x(k  -  2)  4-  w22  wj  -v(A  -  1) 

—  w2  w4  j>(A'  —  1)  —  Wj  y(k  —  2)  4-  Wj  Wj  ;c(A  —  1)  —  Wj  w3  j(A  —  1) 

4-  w4  *(A  -  2)  +  w4  w2  x(k  —  1)  —  w4  w4  >-(A  -  1)  -  w3  y(k  -  2) 

+  w3  Wj  .v(A  —  1)  —  w3  w3  j(A  —  1)  4-  w4  x(A  —  1)  —  w3  y(A  —  1) 

Grouping  terms  we  get, 

y(k)  =  x(k)  4-  (w2  +  w2  w2  4-  w1  w1  4-  w2  w4  4-  vv^  w3  4-  w4)  x(A  —  1) 
4-  (w2  +  w4)  jc(/c  -  2) 

/     1    .        12,        12,        12,        12,        K     /,  1X 

—  (n-;  4-  W4  w2  4-  w3  u-]  +  w4  w4  +  w3  w3  +  w3)>'(A  —  1) 
-(vl.-24-vv2)j(A-2) 


(E-l) 


(£-8) 


From  the  gradient  estimator  and  coefficient  update  equations  we  know  that  the  follow- 
ing relationships  among  lattice  coefficients  are  true 

112  2  112  2  fv       a\ 

Wj  =  w3      Wj  =  vv3      w2  =  w4      w2  =  vv4  {t  —  V) 

Using  these  equalities  in  equation  (E-8),  we  obtain  the  final  expression  for  the  difference 

equation. 


y(k)  =  x{k)  4-  2(w2'  4-  wj  wf  4-  w2  w2)  x(k  -  \)  +  2  w\  x{k  -  2) 
—  2(vvj  4-  w2  w2  +  Wj  W\)y(k  —  1)  —  2  Wj  ^'(A  —  2) 


(£-10) 
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